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SUMMARY 

Solutions  were  obtained  for  hypersonic  flow  over  sharp  cones  at 
high  angle  of  attack  which  include  the  viscous  effects  present  in 
experiment.  These  solutions  were  obtained  by  integrating  the  Navier- 
Stokes  equations  subject  to  a conical  symmetry  assumption.  The  integration 
technique  used  was  MacCormack's  method  with  boundary  conditions  chosen 
to  match  available  experiments.  The  experimental  conditions  were: 

(1)  Those  of  Tracy  for  a 10°  half  angle  cone  at  angles  of  attack  of 
0°,  8°,  12°,  20°,  and  24°  in  M = 7.95  flow  at  Re  = 4.2  x 10^.  (2)  Those 
of  Stetson  for  a 5.6°  half  angle  cone  at  10°  angle  of  attack  in  M = 14.2 
flow  at  Re  = 7.9  x 10^.  (3)  Those  of  McElderry  for  a 6°  half  angle 

cone  at  12°  angle  of  attack  in  M = 6.05  flow  at  Re  = 1.5  x 10^.  A 
physically  based  technique  (normal  stress  damping)  was  demonstrated  for 
controlling  starting  transients  and  for  reducing  or  eliminating  numerical 
oscillations  occurring  at  shock  discontinuities  during  the  integration. 

The  general  features  which  appeared  in  experiment  were  shown  to 
appear  in  the  results  of  the  integration,  including  the  proper  behavior, 
in  laminar  flow,  of  the  viscous  layer  and  the  vortical  singularity.  The 
results  agreed  quite  well  with  the  experimental  data  of  Tracy  (hyper- 
sonic similarity  parameter  x -1)  except  for  a discrepancy  in  pitot 
pressure  in  the  viscous  layer  evident  in  a small  region  near  the  leeward 
centerline  of  the  cone.  The  agreement  with  the  experimental  data  (x  -3) 
of  Stetson  was  less  adequate.  Surface  pressure  agreement  in  this  instance 
was  quite  reasonable.  However,  a somewhat  thin  lee  side  viscous  layer 
resulted  in  a calculated  bow  shock  wave  position  27%  closer  to  the  cone 
surface  at  the  lee  centerline.  It  was  concluded  that  the  lee  side 
viscous  layer  discrepancies  at  both  experimental  conditions  were  primarily 
due  to  lack  of  any  mechanism  in  the  present  technique  to  model  the  non- 
conical  flow  near  the  nose  of  the  cone.  A solution  obtained  just  up- 
stream of  the  beginning  of  boundary  layer  transition  at  the  experimental 
conditions  of  McElderry  agreed  well  with  experiment  when  conically 
projected  far  into  the  turbulent  regine.  The  adequacy  of  the  conical 
symmetry  assumption  is  therefore  indicated  for  the  turbulent  regime  on 
conical  bodies. 
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In  summary,  the  results  show  good  agreement  with  experiment  for 
values  of  the  hypersonic  similarity  parameter  (y  £ 1.0)  and  less  ade- 
quate agreement  at  higher  y-  Angle  of  attack  limitations  encountered 
in  previous  inviscid  cone  flow  calculations  were  not  encountered  in  the 
present  study. 
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SECTION  I 
INTRODUCTION 

Supersonic  flow  about  cones  has  been  a subject  of  interest  to 
aerodynamicists  for  many  years.  Cones  comprised  one  of  the  earliest 
attempts  to  streamline  such  items  as  fuel  tanks,  weaponry,  and  fuse- 
lages, with  locally  conical  nose  shapes  still  in  use  on  today's  modern 
high-speed  aircraft.  The  missile  era  has  also  seen  the  use  of  cones 
as  a primary  or  partial  shape  for  reentry  vehicles  in  current  ICBMs. 

This  study  is  directed  toward  the  latter  application  and  is  in  response 
to  the  need  for  more  accurate  aerodynamic  prediction  techniques  in 
this  area. 

The  two  primary  elements  which  dictate  the  need  for  increased 
prediction  accuracy  are  the  current  high  cost  of  wind  tunnel  testing 
and  the  difficulty  of  simulating  the  reentry  vehicle  flight  reaime. 

The  advent  of  the  high-speed,  large-capacity  computer  in  recent  years 
has  opened  up  a means  to  reduce  cost  both  through  reducina  unit  compu- 
tational cost  and  through  reducing  the  amount  of  experimental  verifi- 
cation required  during  the  systems  acquisition  process.  The  remainina 
task  is  the  research  necessary  to  develop  the  required  numerical  tech- 
niques and  the  verification  of  these  techniques  through  comparison  of 
numerical  results  with  those  obtained  experimentally. 

Particular  solutions  to  the  complete  governing  equations  of  fluid 
flow  are  available  for  very  few  sets  of  boundary  and  initial  condition... 
It  has  therefore  been  necessary  in  nearly  all  cases  to  reduce  the 
complexity  and  scope  of  the  equations  throuah  neolect  of  the  viscous 
effects  and/or  by  reduction  of  the  number  of  dimensions  to  be  considered 
in  the  problem.  Prior  to  1973,  most  attempts  to  solve  the  flow  about 
conical  bodies  used  inviscid  forms  of  the  governinn  equations.  The  tv-n 
primary  difficulties  that  have  been  encountered  with  these  solutions 
are  the  failure  of  the  inviscid  equations  to  model  properly  the  lee  side 
flow  in  which  viscous  effects  predominate  and  the  occurrence  of 
numerical  instabilities  which  are  encountered  at  high  anale  of  attack. 
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If  the  goal  of  reducing  the  amount  of  experimental  verification  is  to  be 
attained,  then  improved  numerical  techniques  must  be  developed  to  over- 
come these  difficulties.  The  present  study  will  present  and  compare 
with  experiment  a technique  which  resolved  both  of  these  difficulties. 

This  report  contains  first  a brief  description  of  the  developments 
in  the  solution  of  flow  over  conical  bodies.  The  solution  approach  used 
in  the  present  study  is  then  presented  in  four  parts.  First,  a discus- 
sion of  the  general  features  of  conical  flows  is  given  with  particular 
emphasis  on  the  circular  sharp  cone  at  angle  of  attack.  The  governing 
equations  are  then  presented  and  changes  due  to  the  conical  symmetry 
assumption  are  discussed.  The  application  of  MacCormack's  numerical 
integration  scheme  to  the  equations  is  given  in  the  third  part.  Part 
four  describes  the  details  of  the  finite  difference  mesh  and  the  manner 
in  which  the  integration  was  accomplished.  The  results  of  the  numerical 
integration  are  then  compared  with  experiment  for  selected  conditions 
for  which  previously  published  experimental  data  were  available. 
Conclusions  drawn  from  these  comparisons  are  presented  to  complete  the 
study. 
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SECTION  II 
BACKGROUND 

The  history  of  the  solution  of  supersonic  cone  flows  has  been  greatly 
affected  by  and  dependent  on  the  development  of  high-speed  digital  com- 
puters. It  is  generally  agreed  that  the  Navier-Stokes  equations  in  their 
complete  form  provide  sufficient  knowledge  of  the  physics  of  the  flow  for 
nonreacting  gases.  However,  the  scope  and  complexity  of  these  partial 
differential  equations  required  great  simplification  of  the  set  in  order 
to  obtain  solutions  prior  to  the  development  of  the  digital  computer. 

These  simplifications  led  to  poor  agreement  of  the  solution  with  experi- 
ment except  for  restricted  conditions.  Although  the  restrictions  have 
been  considerably  relaxed  in  more  recent  studies,  a full  three-dimensional 
solution  of  the  complete  governing  equations  is  still  not  practical  for 
general  configurations.  This  section  will  provide  a brief  description 
of  the  more  important  supersonic  cone  flow  solutions  published  prior  to 
the  present  study. 

The  history  of  solution  techniques  for  supersonic  cone  flow  began 
in  1933  with  the  publication  by  Taylor  and  Maccoll  (Reference  1)  of  an 
ordinary  differential  equation  solution  for  axisymmetric  inviscid  flow 
about  a sharp  cone.  The  next  significant  event  occurred  in  1947  with 
publication  of  tabulated  Taylor-Maccol 1 solutions  by  Kopal  and  Asso- 
ciates (Reference  2).  These  tabulations  were  extended  to  include  angle 
of  attack  in  References  3 and  4.  Ferri  (Reference  5)  obtained  first- 
order  solutions  for  cones  at  small  angle  of  attack  in  which  he  noted 
the  properties  of  the  vortical  layer  about  the  cone  and  showed  the 
existence  of  a vortical  singularity  in  the  lee  side  flow.  With  this 
publication  by  Ferri,  all  important  features  (except  for  the  lee  side 
imbedded  shock  waves)  of  the  inviscid  conical  flow  fields  were  known 
and  defined.  Since  Ferri 's  solution  was  demonstrated  only  for  «/0(._  1-6, 
it  remained  then  to  increase  the  angle  of  attack  for  which  a solution 
could  be  obtained  and  to  account  for  imbedded  shock  waves  and  viscous 
effects  present  in  experiment. 
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Improvement  of  the  digital  computer  in  the  middle  1960's  allowed 
the  solution  of  the  nonlinear  Euler  equations  in  more  complete  form. 
Babenko,  et  al.  (Reference  6)  obtained  smooth-body  solutions  for  the 
inviscid  Euler  equations  in  supersonic  flow  without  further  simp! ification. 
This  work  is  an  early  example  of  the  use  of  finite  difference  techniques 
to  integrate  the  three-dimensional  flow  equations  in  their  nonlinear 
form.  Another  example  of  this  approach  is  the  work  of  Bohachevsky  and 
Rubin  (Reference  7). 

The  finite  difference  scheme  used  in  the  present  work  can  be  traced 
back  to  that  of  Lax  and  Wendroff  (Reference  8).  This  second-order  scheme, 
when  applied  to  the  flow  equations  case  in  conservation  law  form  allows 
shock  discontinuities  to  form  without  special  provisions  during  the 
integration  process.  HacCormack's  (Reference  9)  variant  of  the  Lax- 
Wendroff  scheme  as  applied  by  Kutler  (Reference  10)  extended  the  angle  of 
attack  range  for  which  inviscui  sharp  cone  flows  could  be  obtained. 

Further  applications  (References  11-13)  of  this  technique  to  inviscid 
flows  were  made  by  Kutler,  et  al . to  allow  for  real  gas  effects,  varying 
geometries,  and  passing  through  shock  fronts.  However,  the  inaccuracies 
due  to  neglect  of  viscous  effects  were  not  effectively  dealt  with  until 
the  work  of  Lin  and  Rubin  (References  14,15)  and  Lubard  and  Helliwell 
(Reference  16).  Lin  and  Rubin  solved  a boundary  layer  equation  set 
modified  to  account  for  centrifugal  force  and  cross  flow  diffusion  in  the 
weak  interaction  region  and  a parabolic  set  in  the  hypersonic  tip  region 
to  obtain  solutions  including  viscous  effects  up  to  twice  the  cone  half 
angle.  However,  the  former  approach  requires  input  of  surface  pressure 
from  experiment  or  other  source.  Lubard  and  Helliwell  utilized  an  im- 
plicit finite  difference  technique  similar  to  that  of  Lin  and  Rubin  to 
solve  the  flow  eq'jations  subject  to  a parabolic  assumption  for  the 
streamwise  shear  stress  terms.  This  space  marching  technique  required  an 
accurate  initial  surface  for  the  integration  to  proceed.  The  lack  of 
such  a surface  required  that  the  integration  be  started  at  zero  angle  of 
attack  and  that  the  angle  of  attack  be  gradually  increased  until  the 
desired  value  is  reached.  The  solution  must  then  be  carried  far  enough 
downstream  for  relaxation  of  the  effects  of  this  procedure  to  occur. 

The  comparison  with  experiment  for  this  technique  at  a given  physical 
point  on  the  cone  surface  is  therefore  uncertain. 
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As  noted  in  this  brief  survey,  solutions  of  viscous  flow  over  cones 
have  only  recently  been  attempted  (1973).  The  present  study  seeks  to 
illustrate  a technique  which  includes  viscous  effects  in  cone  flow  cal- 
culations and  removes  limitations  present  in  the  inviscid  calculations. 
The  details  of  this  technique  are  given  in  the  following  section. 
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SECTION  III 
APPROACH 

This  section  presents  the  details  of  the  solution  approach  used  in 
the  present  study.  A brief  description  of  the  general  features  of  coni- 
cal flow  is  presented  in  part  one.  This  is  followed  by  a discussion  in 
part  two  of  the  governing  equations  of  fluid  flow  and  the  changes  in 
these  equations  which  result  from  the  conical  symmetry  assumption.  Part 
three  describes  the  application  of  MacCormack's  numerical  integration 
scheme  to  the  equation  set  obtained  in  part  two.  Part  four  then  gives 
the  details  of  the  solution  procedure. 

1.  DESCRIPTION  OF  CONICAL  FLOW 

As  noted  in  the  introduction  to  the  present  study,  supersonic  flow 
about  conical  bodies  has  been  of  considerable  interest  to  aerodynamicists. 
Supersonic  inviscid  conical  flows  appear  at  first  inspection  to  be  fully 
three-dimensional.  However,  they  are  in  reality  only  two-dimensional 
with  no  gradients  in  any  quantity  occurring  along  the  third  dimension. 

A more  detailed  description  of  these  flows  is  presented  here  to  promote 
understanding  of  their  unique  properties. 

Truly  conical  supersonic  flows  are  steady  state  solutions  to  the 
Euler  equations  (Navier-Stokes  equations  with  viscous  terms  neglected) 
for  a particular  set  of  boundary  conditions.  These  conditions  are  a 
uniform  supersonic  outer  flow  and  a body  generated  by  rays  passing 
through  a common  vertex.  To  ensure  that  the  flow  is  completely  conical, 
the  additional  condition  of  a bow  shock  wave  attached  to  the  vertex  of 
the  body  must  be  imposed.  Examples  of  bodies  which  generate  these  flows 
are:  (1)  circular  cones,  (2)  elliptical  cones,  (3)  delta  wings, 

(4)  delta  wing/conical  body  combinations,  and  (5)  axial  corners.  Many 
variations/combinations  of  these  bodies  are  possible  which  produce 
conical  flows.  Also,  conical  bodies  which  appear  as  nose  shapes  in  more 
general  configurations  will  produce  locally  conical  flows.  The  char- 
acteristics of  these  flows  are  examined  below. 
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The  dominant  feature  of  conical  flow  is  that  all  fluid  quantities 
are  constant  along  rays  passing  through  the  vertex  of  the  conical  body. 

By  noting  that  these  rays  are  equivalent  to  the  spherical  radius  r with 
origin  at  the  cone  vertex,  conical  flow  can  be  mathematically  described 
by  stating  that  derivatives  of  all  fluid  quantities  with  respect  to  r are 
identically  zero.  The  physical  implications  of  this  statement  are 
illustrated  in  Figure  1.  This  figure  shows  a plane  cut  of  an  axisymmetric 
supersonic  cone  flow.  A conical  flow  exists  when  the  flow  quantities  at 
point  A (along  streamline  are  identical  to  those  at  point  B (along 
streamline  fg).  The  flow  at  both  point  A and  point  B can  be  completely 
described  in  terms  of  the  single  angle  9 for  this  two-dimensional  cut. 

What  is  apparently  a two-dimensional  flow  is  then  in  reality  only  one- 
dimensional (sometimes  referred  to  as  1-1/2  dimensional,  since  two 
velocity  components  are  still  present). 

Although  somewhat  more  difficult  to  visualize,  conical  bodies  at 
angle  of  attack  in  supersonic  flow  also  produce  conical  flow  fields.  By 
analogy  with  the  above  discussion,  these  flows  can  also  be  described 
completely  in  terms  of  the  independent  variable  0 and  a circumferential 
angle  4>  (these  angles  are  illustrated  in  Figure  2).  All  natural  features 
of  the  flow;  such  as  imbedded  shocks,  bow  shocks,  slip  lines,  etc.;  form 
surfaces  composed  of  rays  passing  through  the  conical  vertex.  An  example 
sketch  of  this  flow  is  shown  in  Figure  3.  The  streamlines  shown  in 
Figure  3 show  the  manner  in  which  the  flow  crosses  the  bow  shock,  expands 
around  the  body,  and  then  is  turned  back  along  the  cone  by  the  imbedded 
shock.  It  is  instructive  to  view  flow  features  on  a spherical  surface 
through  the  field. 

Figure  4 illustrates  the  approximate  cross  flow  streamline  pattern 
on  the  spherical  surface  and  shows  relative  location  of  the  shock  waves. 
The  point  at  which  the  cross  flow  streamlines  converge  is  the  vortical 
singularity  first  noted  by  Ferri  (Reference  5).  The  streamlines  all 
have  different  values  of  entropy  which  implies  that  density  and  entropy 
are  discontinuous  at  the  singularity  with  pressure  beina  continuous. 
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This  singularity  rests  on  the  cone  surface  (in  fact,  Melnik  (Reference 
17)  has  shown  the  possibility  of  two  singularities  on  the  surface)  at  low 
angles  of  attack  and  lifts  off  as  angle  of  attack  increases. 


In  summary,  conical  flows  can  be  completely  described  in  a spherical 
coordinate  system  by  only  two  independent  coordinates  (0,  <j)).  All 
natural  features  of  the  flow  are  also  conical  surfaces.  These  facts 
tend  to  render  a very  complex  flow  more  amenable  to  solution  by 
presently  available  techniques. 


2.  THE  FLUID  FLOW  EQUATIONS 

In  this  section,  the  basic  equations  of  fluid  flow  are  presented 
in  complete  form  and  the  conical  symmetry  assumption  is  applied.  The 
lesulting  equations  allow  integration  to  take  place  on  a single  spherical 
surface.  The  non-dimensional ization  of  the  equations  and  the  physical 
.I'eaning  of  the  assumptions  are  discussed. 

The  Navier-Stokes  equations  which  describe  flow  of  a perfect  gas, 
ran  be  written  in  conservation  law  form  (Reference  18)  for  a spherical 
coordinate  system  as  follows: 
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The  stress  tensor  terms  and  o^-^.  will  be  defined  in  detail  later. 

Solution  of  the  equation  set  in  the  form  shewn  above  would  require 
very  large  amounts  of  machine  time  even  for  simple  aerodynamic  config- 
urations. This  has  led  researchers  in  the  past  to  simplify  the  equation 
set  through  neglect  of  the  stress  terms  (inviscid  flow)  and/or  by 
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reducing  the  number  of  dimensions  to  be  considered  in  the  problem.  One 
class  of  flows  which  has  been  extensively  examined,  (References  10,19,20) 
through  neglect  of  the  viscous  terms,  is  that  of  inviscid  conical  flow. 


A conical  flow,  as  noted  in  the  previous  section,  can  be  described 
as  an  inviscid  steady  flow  in  which  all  flow  quantities  are  constant  along 
rays  passing  through  the  vertex  of  the  conical  body.  If  a body  fixed 
spherical  coordinate  system  (Figure  2)  centered  at  the  vertex  of  the 
conical  body  is  used  to  describe  the  flow,  then  all  spherical  surfaces 
must  have  the  same  vector  and  scalar  values  of  the  flow  quantities  for  a 
given  (0,  <{))  point  on  each  surface.  Therefore,  all  derivatives  of  flow 
quantities  with  respect  to  the  spherical  radius  (r)  of  these  spherical 
surfaces  from  the  origin  must  be  zero.  This  has  the  effect  of  reducing 
the  number  of  independent  variables  in  the  problem  by  one. 


Examination  of  experimental  studies  of  supersonic  flow  over  conical 

bodies  (References  21-23)  reveals  that  these  ■•'lows  exhibit  approximate 

conical  behavior  downstream  from  the  nose  region  even  though  relatively 

large  viscous  regions  exist.  Cross  (Reference  21)  determined  that  the 

viscous  layer  growth  on  the  lee  side  of  a delta  wing  in  supersonic  flow 

was  essentially  conical.  The  oil  flow  separation  traces  for  sharp  cones 

in  the  experimental  study  of  Stetson  (Reference  22)  are  approximately 

straight  (but  the  conical  vertex  of  these  traces  is  displaced  downstream 

by  nose  effects).  Therefore,  in  concert  with  an  idea  first  broached  by 

Anderson  (Reference  24)  for  axial  corner  flow,  the  conically  symmetric 

flow  approximation  (-^:^  HO)  is  applied  to  all  terms  in  Equation  1. 
o r 


The  resulting  equation  set  is: 
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where; 

(a)  d,  f,  and  S are  unchanged  except  in  the  definition  of  the 
stress  terms. 
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These  equations  are  non-dimensional ized  as  follows,  with  dimensional 
quantities  denoted  by  primes: 


u = u_ 


V = _y_ 
V„, 


w = _w_ 
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p = p = t = r 


T = II 


n V P4.  r 
Re  = max  t°° 


e = e' 
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V = f2y  R t;  1 
kfr  9 


Equation  3 indicates  a rotational  change  from  this  point  onward. 
Hereafter,  units  will  be  specified  where  dimensional  quantities  are 
used.  The  quantities  used  for  non-dimensional ization  are  due  to 
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Kutler  (Reference  10).  it  is  significant  that  when  time  is  non-dimension- 
al ized  by  the  parameter  the  r dependence  in  Equation  2 is 

contained  in  the  Reynolds  number.  The  net  effect  is  that  the  calculation 
is  carried  out  at  a single  (6,  (}>)  spherical  surface  with  the  distance  of 
this  surface  from  the  cone  apex  determined  by  the  Reynolds  number. 
Therefore,  all  spherical  radius  scaling  is  now  contained  in  the  Reynolds 
number  alone. 


as  follows; 
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The  heat  transfer  terms  are  defined  as  follows: 


^9 


2RePr  36 


2RePr  3;}) 


where 


Pr  = — p. 

K 


(5) 


At  this  point,  it  is  prudent  to  examine  further  the  physical 
implications  of  the  approximations  applied  to  the  above  equation  set. 
It  should  first  be  noted  that  the  only  approximations  now  inherent  to 
the  equations  are  that  a perfect  gas  is  required  and  that  conically 
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symmetric  flow  is  assumed  at  all  points  in  the  field.  Figure  5 depicts  a 
cross  section  through  an  axi symmetric  cone  flow  with  a somewhat  exaggerated 
boundary  layer  thickness  shown  by  the  dashed  line.  The  thesis  is  that  at 
a given  calculation  surface,  the  viscous  layer  thickness  will  be  properly 
scaled  by  the  Reynolds  number  based  on  the  spherical  radius  to  that 
surface.  The  only  error  then  would  be  in  the  local  gradients  of  the  flow 
quantities  and  would  be  of  the  order  of  the  angle  between  the  spherical 
radius  and  the  edge  of  the  boundary  layer  at  the  calculation  surface. 

The  same  statement  could  be  made  concerning  a calculation  at  any  point 
along  the  cone  downstream  from  the  non-conical  nose  region.  Therefore, 
since  the  growth  of  the  viscous  layer  thickness  will  not  in  general  be 
linear,  the  assumption  in  reality  only  implies  locally  conical  flow. 

In  order  to  test  the  above  hypothesis,  calculations  were  carried 
out  using  the  above  equation  set  and  the  method  set  forth  in  the 
following  sections  applied  to  axisymmetric  cone  flow.  The  cone  surface 
pressure  resulting  from  these  calculations  is  shown  in  Figure  6 and 
compared  with  hypersonic  weak  interaction  theory  for  axisymmetric  cone 
flow  (Reference  29).  The  agreement  shown  in  Figure  5 between  the  cal- 
culations and  the  weak  interaction  theory  is  quite  good  with  divergence 
noted  (as  expected)  toward  the  nose  of  the  cone.  In  order  to  fully 
verify  the  hypothesis,  two-dimensional  calculations  compared  with 
appropriate  experimental  data  are  required.  The  remainder  of  this  report 
will  present  the  results  of  calculations  carried  out  at  selected  con- 
ditions for  which  experimental  data  are  available.  Comparisons  are 
made  and  conclusions  are  drawn  concerning  the  adequacy  of  the  hypothesis. 

3.  THE  NUMERICAL  INTEGRATION  SCHEME 

This  section  discusses  the  choice  of  MacCormack's  finite  difference 
scheme  as  the  integration  method  for  this  problem.  The  scheme  is  il- 
lustratively applied  to  the  one-dimensional  wave  equation  as  an  example. 

The  equations  to  be  integrated  are  cast  in  MacCormack's  predictor- 
corrector  form,  and  the  finite  difference  representation  of  the 
derivatives  is  shown. 
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The  utility  of  MacCormack's  (Reference  9)  finite  difference  scheme 
for  the  solution  of  in/iscid  supersonic  flows  has  been  demonstrated  in 
numerous  studies  (References  10-13).  More  recent  studies  have  also  shown 
it  to  give  excellent  results  for  two-dimensional  viscous  flows  with  sepa- 
ration (References  25-27).  The  scheme  can  be  described  as  a second-order 
accurate  predictor-corrector  sequence  for  the  integration  of  partial 
differential  equations.  MacCormack's  scheme  is  a variant  of  the  Lax- 
Wendroff  scheme  (Reference  10)  and  can  be  shown  (as  applied  to  the 
linear  wave  equation)  to  reduce  to  Lax-Wendroff  form  when  the  predictor  is 
substituted  into  the  corrector.  Since  second-order  difference  schemes 
generally  give  better  results  with  fewer  mesh  points  than  first-order 
schemes,  MacCormack's  proven  differencing  scheme  is  chosen  for  the 
present  study. 


If  an  Euler  predictor  followed  by  a modified  Euler  corrector  is 


applied  to  the 

1-D  wave  equation 

3U  ^ 3U  _ „ 
at  ^3x  " 

Predictor: 

»i+i  = - 

CAt  (U^). 

(6) 

Corrector: 

+ - CAt  (&^) 

(7) 

MacCormack's  scheme  results  when  forward  first  spatial  differences  are 
used  in  the  predictor  and  backward  first  differences  in  the  corrector. 
The  predictor  corrector  sequence  is  then  (with  i representing  time 
location  and  j spatial  location): 


Predictor : 


Corrector: 
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The  two  equations  collapse,  when  the  appropriate  substitutions  are  made, 
to  the  following  Lax-Wendroff  form: 
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The  leading  error  terms  which  are  truncated  through  use  of  this  finite 
difference  equation  representation  are  of  the  form  (with  the  derivatives 
evaluated  at  an  undetermined  point  within  the  range  of  interest) 

^ cAtAx^  ,,  cAt^Ax^  „ , At^  „ 

~6  ^x^  • —24  ^x^  ^ ^t’ 

which  illustrates  that  when  At  and  Ax  are  of  the  same  order,  the 
estimated  truncation  error  must  be  of  0(At  ) and  the  method  can  be 
said  to  be  second-order  accurate.  The  stability  of  MacCormack's  scheme 
is  examined  in  Appendix  C. 

When  the  scheme  is  applied  to  Equation  2,  the  resulting  MacCormack's 
predictor-corrector  steps  are  as  follows  (the  » used  previously  to 
indicate  vector  quantities  is  omitted  here  for  clarity): 


Predictor: 


Corrector : 


where  n.  indicates  that  the  flow  quantities  are  evaluated  at  the  predictor 
level  and  o]  . implies  D(iAt,  jA0,  kA<{)).  The  presence  of  the  H'  matrix 
in  these  equations  indicates  that  they  are  in  so  called  "weak"  conservation 
form  (Reference  28).  Numerical  results  in  the  present  study  reveal  that 
no  significant  error  in  total  temperature  occurs  across  shock  transitions 
through  use  of  this  form  of  the  equations  (implying  that  energy  must  have 
been  properly  conserved). 


Since  Equation  2 contains  derivatives  in  the  stress  terms  which 
remain  as  an  integral  part  of  the  matrix  terms,  a value  of  the  derivatives 
must  be  obtained  at  each  of  the  mesh  points  used  in  the  primary  differencing. 
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inese  values  must  be  obtained  in  a manner  which  will  maintain  both 
consistency  and  accuracy  of  the  differencing  equation.  This  was  done 
in  the  present  study  by  evaluating  the  interior  derivative  terms  with 
respect  to  the  same  independent  variable  using  first-order  difference 
approximations  of  opposite  sense  to  the  exterior  differences.  An 
illustration  follows  for  the  predictor. 


Consider  the  vector  which  consists  of  functions  of  the  flow 
quantities  and  their  derivatives.  The  finite  difference  representation 

sr 

30 


Of 


1 s 


3£ 

30 


A6 


(13) 


Assume  that  the  vector  1^  contains  the  derivative  values  of  which  are 

therefore  required  at  mesh  points  j,k  and  j+l,k  in  the  difference 
representation.  These  values  are  obtained  by  ^*^3+1  k * k^ 

at  j-tl,k  and  - Izr  (o,-  , - u.  , .)  at  j,k.  The  resulting  finite 

do  At)  J ^ K j * ‘ 9 K 

difference  representation  for  is  the  standard  second-order  accurate 
difference  centered  at  j. 

When  a cross  derivative  term  is  involved  (say  a term  contained 
^ <•(}) 

in  T— ) , second  order  centered  first  differences  are  used.  With  -- 
30  3u  1 

represented  as  above,  at  j+l,k  and 

(Uj  - Uj  |^_.|  ) at  j,k.  This  results  in  the  standard  second- 

^ u 

order  accurate  second  difference  for  the  difference  centered  at  j+1/2. 

d(fi30 

This  procedure  is  also  used  in  the  corrector  step  with  all  differencing 
reversed  in  the  0 direction.  The  resulting  differences,  when  combined 
with  the  predictor  into  the  Lax-Wendroff  form,  are  correct  and  centered 
at  j. 


This  essentially  completes  the  application  of  MacCormack's  scheme 
to  Equation  2.  It  should  be  apparent  that  throughout  the  analysis  of 
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MacCormack ' s scheme,  careful  prescription  of  the  first-order  differences 
in  the  individual  predictor  and  corrector  results  in  a second-order 
accurate  scheme  for  the  entire  cycle.  The  details  of  the  integration 
procedure  are  described  in  the  following  section. 

4.  THE  SOLUTION  PROCEDURE 

In  this  section,  the  details  of  the  integration  procedure  are 
discussed.  The  finite  difference  mesh  is  first  described,  followed  by 
the  initial  i/:ation  used  for  the  initial  value  problem  and  a brief 
description  and  discussion  of  the  physical  boundary  conditions.  Finally, 
the  convergence  criteria  for  determining  when  the  sought-for  steady 
state  solution  has  been  reached  is  described. 

As  stated  earlier  in  this  section  the  calculation  takes  place  on  a 
spherical  surface  located  at  a distance  r from  the  vertex  of  the  cone. 

The  distance  r is  determined  by  the  distance  along  the  cone  at  which 
correlation  with  experiment  is  desired.  It  then  appears  as  the  character- 
istic length  in  the  Reynolds  number  and  scales  the  viscous  effects  at  the 
calculation  surface.  The  integration  procedure,  as  with  all  finite 
difference  techniques,  requires  values  of  the  Fluid  flow  quantities  to  be 
known  at  discrete  points  labeled  as  mesh  points.  For  the  present  study, 
these  points  are  arranged  on  the  (0,  f)  calculation  surface  as  shown  in 
Figure  2.  The  mesh  points  are  evenly  spaced  in  the  0 and  $ directions 
with  A6  f A((i.  It  should  be  noted  that  considerable  advantage  for 
viscous  calculations  can  sometimes  be  obtained  by  varying  mesh  point 
spacing  to  cluster  points  near  the  surface.  However,  only  uniform 
size  mesh  is  considered  here  since  the  present  study  is  primarily 
concerned  with  proof  of  an  untried  assumption  for  the  flow  equations. 
Therefore,  the  additional  complexity  involved  in  obtaining  the  uneven 
mesh  spacing  was  not  desirable  and  was  not  made  a part  of  this  study. 

The  mesh  spacing  used  was  designed  to  give  entirely  adequate  results 
on  the  lee  side  of  the  cone  in  the  large  viscous  regions  present  there. 

The  finite  difference  mesh  is  initialized  to  free  stream  values 
of  the  flow  quantities  except  for  surface  points  at  which  the  three 
velocity  components  are  set  identically  to  zero.  This  is  physically 
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equivalent  to  bringing  the  cone  from  rest  to  the  free  stream  Mach 
number  instantaneously.  As  might  be  expected,  the  initial  transients 
associated  with  this  procedure  are  quite  large  and  require  numerical 
damping  at  the  higher  angles  of  attack.  However,  this  is  not  considered 
detrimental  since  the  damping  is  required  continuously  for  most  of  the 
solutions  obtained  for  reasons  which  will  be  discussed  later.  Also,  a 
comparison  between  final  results  obtained  with  the  impulsive  start  and 
with  a solution  obtained  by  changing  slightly  the  boundary  condition 
representation  (thereby  perturbing  the  flow)  of  a converged  solution 
revealed  no  essential  difference  between  the  two  solutions. 

The  boundary  conditions,  which  are  maintained  throughout  the 
calculation,  are: 

(1)  All  velocities  identically  zero  at  the  surface. 

(2)  Surface  temperature  constant  and  equal  to  the  value  for  the 
experimental  case. 

(3)  Free  stream  values  of  all  flow  quantities  are  maintained  at 
the  outer  edge  of  the  mesh.  Since  the  bow  shock  wave  was  very  close 
to  the  cone  surface  on  the  windward  side  for  the  high  < cases,  an 
ellipse  was  used  to  provide  a dummy  outer  edge  of  the  mesh  which  would 
be  close  to  the  bow  shock  and  therefore  save  computer  time. 

(4)  Lateral  symmetry  of  the  flow  quantities  is  maintained  across 
the  t = 0°  and  = 180®  lines. 

The  values  of  pressure  and  density  at  the  surface  are  determined  during 
the  integration.  A discussion  of  this  matter,  plus  a more  complete 
description  of  the  boundary  conditions,  appears  in  Appendix  A. 

In  this  study,  the  predictor  and  corrector  are  swept  through  the 
((fi,  O)  mesh  in  turn,  beginning  at  the  .J)  = 0°,  0 = 6^  point.  A complete 
traverse  of  the  predictor  and  corrector  constitutes  one  time  step. 

Step  by  step,  the  procedure  is  as  follows: 

(l)  Tne  predictor  is  swept  through  the  mesh  from  9 = 0^  to  0 = 0^ 
and  from  $ = 0°  to  180®. 
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(2)  The  l5  matrix  is  decoded  to  obtain  updated  values  of  the  flow 
quantities. 

(3)  The  corrector  is  swept  through  the  mesh,  using  the  updated 
flow  quantities  where  appropriate. 

(4)  The  ^ matrix  is  decoded  to  give  values  for  the  flow  quantities 
at  the  end  of  the  time  step. 

(5)  Boundary  conditions  are  applied  where  necessary  during  the 
integration. 

As  implied  in  the  stated  outer  boundary  condition,  no  special 
provisions  are  included  to  handle  the  bow  shock  wave  or  any  imbedded 
shock  waves  which  may  arise.  These  shock  waves  are  "captured"  through 
use  of  the  conservation  form  of  the  governing  equations  (References  10,28). 
This  procedure  is  considered  adequate  through  use  of  numerical  damping 
for  reducing  spurious  oscillations  around  the  shock  waves  (Appendix  B). 

A criteria  for  determining  when  the  steady  state  has  been  reached 
(i.e.  "convergence")  is  considered  an  integral  and  necessary  part  of 
this  study.  Previous  unpublished  experience  by  the  author  with  inviscid 
flows  reveals  that  instabilities  in  the  calculation  can  arise  even 
after  changes  in  the  flow  field  pressure  become  small  as  seen  on  a 
plotter  or  CRT  display  device.  Therefore,  a stringent  convergence 
criteria  is  used  in  the  present  study  in  order  to  insure  that  the 
solution  is  stable  and  that  the  results  obtained  are  repeatable.  The 
convergence  criteria  that  was  used  required  computation  of  differences 
in  pressure,  density,  velocity,  and  energy  between  time  steps  and 
stopped  the  calculation  when  these  differences  reached  the  fifth  signi- 
ficant figure  for  all  points  in  the  flow  field.  At  the  time  step  sizes 
used  for  runs  reported  here,  convergence  of  the  solution  typically 
occurred  at  the  physical  time  required  for  the  outer  flow  to  move  four 
to  five  times  the  length  from  the  nose  to  the  calculation  surface.  By 
reducing  the  number  of  significant  figures  input  to  the  convergence 
routine,  it  was  found  that  the  fifth-place  criteria  gave  run  times 
nearly  twice  as  long  as  would  have  been  required  to  obtain  values  of 


19 


AFFDL-TR-76-139 


surface  pressure  accurate  to  engineering  standards  (Changes  of  less  than 
3%  occurred  in  surface  pressure  during  the  last  50%  of  the  run).  The 
numerical  stability  of  the  present  calculations  is  therefore  assured. 

The  present  study  did  not  make  continuous  use  of  a stability  limit 
criteria  to  determine  allowable  step  size.  Instead,  the  time  step  size 
for  each  computer  run  was  frozen  at  a value  which  would  insure  no 
difficulty  with  normal  stress  damping  (described  in  Appendix  B).  Based 
on  experience  gained  in  the  study,  a stability  criteria  was  postulated 
and  confirmed  by  numerical  experiment.  This  entailed  repeat  of  a 
previous  run  with  time  step  size  automatically  determined  by  applying 
the  criteria  noted  below  at  each  time  step.  This  stability  limit  criteria 
is  described  in  detail  in  Appendix  C.  The  criteria  is  summarized  here: 


At  = m MIN  (At^^^,  At^^^,  At^.^^) 


(14) 


where  m < 1 


At 


cfl-  _[vi 
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A9(sin0A<}) 


iV2 
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pRe  A9^ 
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A9^  (fRe  Pr) 
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In  this  instance,  MIN  implies  the  minimum  value  of  the  quantity 
found  by  searching  all  mesh  points  in  the  flow  field.  Although  m is 
theoretically  limited  to  1.0  (the  CFL  condition),  it  was  found  that  a 
value  of  m = 1.2  could  be  used  once  the  initial  strong  transients  were 
past.  Instability  of  the  solution  occurred  at  m = 1.4. 


The  operation  sequence  used  in  the  computer  code  was  the  following: 


0) 

READ  input  data 

(2) 

Initial ize  mesh 

(3) 

DO  time  steps,  1 

to  input  number 

(4) 

Calculate  local 

Reynolds  number 
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(5)  Calculate  local  damping  coefficients 

(6)  DO  predictor  for  all  but  free  stream  boundary  mesh  points 

(a)  Calculate  coordinate  metrics 

(b)  CALL  boundary  value  subroutine 

(c)  Define  finite  difference  representations 

(d)  CALL  vector  load  subroutine  (creates  5,  and  in 
Equation  2 according  to  the  requirements  of  (c)  above) 

(e)  CALI  integrations  subroutine  (creates  predicted  value 
of  0 at  new  time  level  by  application  of  Equation  11 

(f)  CALL  [5  vector  solver  subroutine  (decodes  flow  quantities 
from  calculated  l5  vector  at  new  time  level) 

(7)  CALL  symmetry  boundary  condition  subroutine 

(8)  DO  corrector  for  all  but  free  stream  boundary  mesh  points 
(This  is  an  exact  duplicate  of  step  (6)  except  for  changes  in  finite 
difference  representations  and  flow  quantity  time  levels.  The  code 
used  the  same  subroutines  for  both  predictor  and  corrector  sequences.) 

(9)  CALL  symmetry  boundary  condition  subroutine 

(10)  Test  for  convergence.  If  not  converged,  return  to  step  (3)  and 
continue  until  convergence  criteria  is  met  or  the  input  number  of  time 
steps  is  completed. 

(11)  Store  results 

(12)  Print  and  plot  results  of  calculation 

(13)  STOP 

This  sequence  of  operations  may  be  used  to  integrate  (by  use  of 
MacCormack's  scheme)  other  three  coordinate  versions  of  the  flow 
equations  by  appropriate  changes  in  steps  (6)  (b),  and  (d).  The 
restrictions  are  that  weak  conservation  form  be  used  and  that  the 
independent  variable  of  integration  be  time  or  a coordinate  direction 
in  which  M > 1 . 
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This  completes  the  description  of  the  approach  used  in  the  present 
study.  The  following  section  describes  and  compares  with  experimental 
. data  the  numerical  results  obtained  through  application  of  this  approach. 


\ 
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SECTION  IV 
NUMERICAL  RESULTS 

This  section  presents  a comparison  between  results  obtained  with  the 
present  technique,  the  numerical  technique  of  Lubard  and  Helliwell,  and 
selected  wind  tunnel  experiments.  The  first  comparison  is  with  the 
experimental  study  of  Tracy  (Reference  30).  This  study  has  been  a 
standard  of  comparison  for  other  techniques  (References  14-16)  and  is 
very  well  documented  in  the  open  literature.  The  second  comparison  is 
with  the  experimental  study  of  Stetson  (Reference  22)  and  the  third 
with  that  of  McElderry  (Reference  31). 

The  validity  of  any  heretofore  untried  assumption  applied  to  the 
governing  equations  can  only  be  determined  through  comparison  of  the 
results  obtained  by  calculation  with  results  of  experiment.  To  test 
the  validity  of  the  present  technique,  experimental  data  were  chosen 
which  had  been  published  and  which  provided  a Mach  number  range  of  six  to 
fourteen.  The  experimental  cases  chosen  were  for  sharp  cones  within 
the  limits  reasonably  obtainable  by  standard  machining  techniques  (0.001 
to  0.002  in.  nose  radius).  Reynolds  numbers,  except  for  that  of  McElderry, 
were  in  the  range  for  which  laminar  flow  could  be  expected.  Transition 
onset  in  McElderry's  data  was  at  approximately  Re  = 2.4  x 10^.  The 
cases  all  have  flow  field  probe  data  available  for  selected  Reynolds 
numbers  and  angles  of  attack.  The  remainder  of  this  section  will  present 
a detailed  comparison  of  the  results  of  the  present  calculation  with  the 
chosen  experimental  data. 

1.  TRACY'S  DATA 

Tracy  (Reference  29)  conducted  an  experimental  study  using  a 10  half 
angle  sharp  (0.002  nose  radius)  cone  in  M = 7.95  flow.  The  wind  tunnel 
working  fluid  was  air  with  total  pressure  and  temperature  of  259.3  psia 
and  1360''R  for  the  runs  of  interest.  The  resulting  free  stream  Re  was 
1.25  X 10^/ft.  Data  were  taken  at  angles  of  attack  of  O'  to  24°  in  ^ 
increments.  Measurements  of  surface  pressure  at  x = 4.0  in.  and  surveys 
of  pitot  pressure  at  x - 3.45  in.  (and  perpendicular  to  the  cone  surface) 
were  selected  for  comparison.  The  conf iqurat ion  of  the  model  and 
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the  same  x station  on  the  cone  surface. 

Calculations  were  performed  at  the  following  conditions; 

1 

M = 7.95  Re  = 4.2  x 10^  (x  = 4.0  in.)  [ 

X 

= 259.3  psia  ; 

00 

T.  = 1360°R  T/T,  =0.41 

t w t 

CO  00 

I a = 0°,  8°,  12°,  20°,  20.5°,  24°  Pr  = 0.72 

The  finite  difference  mesh  contained  48  x 50  points  except  for  the  a = 24°  j 

case  in  which  a 48  x 70  mesh  was  used.  The  circumferential  ((*))  step  size  i 

was  4°  for  all  runs  with  the  0 step  size  varying  depending  on  the  angle  of  i 

attack  and  number  of  mesh  points  used.  All  runs  in  this  study  were  con-  j 

verged  to  the  fifth  significant  digit  of  the  flow  quantities.  Computer 

run  times  for  these  solutions  ranged  from  1.4  hours  (a  = 12°)  to  1.6  hours  \ 

(a  = 24°)  on  the  CDC  CYBER  74  (equivalent  in  CPU  speed  to  a CDC  6600). 

The  number  of  time  steps  required  ranged  from  1860  (a  = 8°)  to  2480  (a  = 12°). 

The  fact  that  the  highest  angle  of  attack  run  did  not  take  the  largest 
number  of  time  steps  is  due  to  variance  in  the  actual  CFL  condition  at  which 
each  calculation  was  made.  More  consistent  results  are  obtained  through 
use  of  the  stability  criteria  set  forth  in  Appendix  C.  In  order  to 
provide  comparison  of  both  flow  field  and  surface  quantities,  comparison 
with  experimental  data  is  shown  for  surface  pressure  and  nitot  pressure 
surveys.  In  addition,  velocity  vector  plots  are  shown  for  the  calculations 
and  pictorial  displays  of  cross  flow  field  features  are  shown  ’!n  order 
to  compare  the  overall  agreement  of  the  calculation  with  experiment. 

(These  data  were  also  reported  in  (Reference  39). 

Figures  7 through  10  compare  the  numerical  results  of  the  calculation 
with  the  surface  pressure  data  of  Tracy.  At  this  point  it  should  be 
noted  that,  unless  otherwise  stated,  the  discrete  points  shown  for 
Tracy's  data  were  obtained  from  continuous  curves  in  Reference  30 
through  use  of  a digitizer.  Figure  7 gives  this  comparison  for  a = 8°. 

The  calculated  value  can  be  seen  to  fall  just  at  the  edge  of  the  symbols 
for  the  entire  circumference  of  the  cone.  This  corresponds  to  a local 
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difference  of  1.67'  at  = 0°,  the  windward  centerline  of  the  cone. 

The  surface  pressure  comparison  for  -i.  = 12°  is  given  in  Figure  8.  The 
agreement  in  this  case  is  somewhat  better  near  the  leeward  centerline 
(({>  = 180°)  but  the  local  difference  has  increased  to  2.5"  at  <})  = 0. 

The  same  type  of  agreement  exists  for  the  a = 20°  (Figure  9)  and  a = 24° 
(Figure  10)  cases  with  the  local  difference  for  a = 24°  becoming  5%  at 
4>  = 0°.  To  provide  a further  comparison,  the  surface  pressure  at  $ = 0° 
obtained  by  Lubard  and  Helliwell  (Reference  16)  for  the  a = 12°  case  is 
identical  to  that  obtained  in  the  present  study.  The  first  possible 
reason  for  the  discrepancy  between  calculation  and  experiment  on  the 
windward  side  is  that  Tracy  noted  an  approximate  0.5°  uncertainty  in 
angle  of  attack  in  the  wind  tunnel  used  for  the  experiment.  In  order 
to  test  this  possibility,  the  a = 20°  run  was  reconverged  at  a = 20.5°. 
The  surface  pressure  at  cj)  = 0°  regained  only  1/2  of  the  4.7%  difference 
between  the  experiment  and  the  a = 20°  solution,  indicating  that  not 
all  of  the  discrepancy  could  be  due  to  the  uncertainty  in  angle  of 
attack.  The  second  possibility  which  should  be  considered  is  the  fact 
that  the  surface  pressure  taps  used  for  the  experiment  were  somewhat 
large  in  relation  to  boundary  layer  thickness  on  the  windward  side. 

In  any  case,  an  error  band  of  at  least  5%  is  present  in  most  experimental 
studies . 

Sample  circumferential  pitot  surveys  are  presented  in  Figures  11 
through  13.  The  angle  in  the  9 direction  for  comparison  with  experiment 
is  defined  by 

9 = 9^  + tan"^  (15; 

where  y is  the  normal  distance  above  the  cone  surface  in  inches.  The 

5 

Reynolds  number  for  the  calculated  results  is  4.2  x 10' . Since  the 
pitot  surveys  were  made  at  constant  height  above  the  cone  surface, 
they  will  in  most  cases  intersect  the  bow  shock  at  some  point  in  the 
circumferential  traverse.  This  appears  in  the  pitot  pressure  plots  as 
a sharp  rise  ■•"rom  the  free  stream  value  (8.7  x 10  "^)  to  the  maximum 
value  for  each  trace.  Probe  effects  cause  the  character  of  the  rise  of 
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experimental  pitot  pressure  to  appear  remarkably  like  that  of  the  rise 
of  pitot  pressure  through  the  captured  shock  waves  obtained  in  the 
present  study,  with  decrease  in  transition  slope  and  loss  of  peak  pressure 
in  the  calculations  caused  by  the  necessary  damping  present.  Note  that 
the  pitot  pressure  transition  centerpoint  (the  shock  location  criteria 
used  by  Tracy)  of  the  captured  shock  in  all  of  the  calculated  results 
is  virtually  the  same  as  would  be  determined  from  the  experimental  pitot 
survey  and  that  the  calculated  pitot  pressure  quickly  returns  to  the 
proper  magnitude  on  the  high  pressure  side  of  the  shock.  For  the  a = 12° 
case  (Figure  11),  the  uppermost  survey  (y  = 0.4  in.)  is  entirely  in  the 
inviscid  flow  region  for  both  calculation  and  experiment.  The  survey 
nearest  the  cone  surface  (v  = 0.05  in.)  is  in  the  viscous  layer  for  more 
than  90°  of  the  circumference  of  the  cone.  Agreement  with  experiment 
for  both  of  these  surveys  is  good.  The  calculated  center  survey  (y  = 0.2 
in.)  departs  from  the  experimental  value  near  the  lee  centerline  (cj)  = 180°). 
This  departure  is  due  to  the  smaller  viscous  layer  present  in  the  calcu- 
lation as  compared  to  experiment. 

The  a = 20°  case  (Figure  12)  has  essentially  the  same  degree  of 
agreement  between  the  calculated  results  and  experiment  with  the  exception 
that  the  departure  near  the  lee  centerline  is  concentrated  in  a region 
of  approximately  6"  (tp  - 174°  - 180°).  Agreement  outside  this  region  is 
generally  good.  Figure  13  illustrates  that  essentially  the  same 
agreement  and  regidVi  of  discrepancy  exist  for  the  a = 24°  case. 

The  pitot  pressure  discrepancy  near  the  lee  centerline  can  be 
evaluated  more  clearly  from  Figure  14,  which  compares  pitot  pressure 
along  the  lee  centerline  for  the  a = 24°  case.  The  discrepancy,  as 
noted  above,  extends  for  approximately  6°  to  either  side  of  the  lee 
centerline  for  this  angle  of  attack.  It  is  attributed  to  a combination 
of  locally  large  pitot  pressure  drop  due  to  finite  model  nose  radius 
and  to  nose  effects  resulting  from  the  non-conical  region  at  the  nose. 

Any  persistent  experimental  nose  effects  would  tend  to  be  concentrated 
along  the  lee  centerline  as  they  are  swept  back  along  the  cone.  This 
phenomenon  can  be  observed  in  the  experimental  study  of  Stetson 
(Reference  22)  through  examination  of  the  oil  flow  pictures  presented 
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therein.  For  instance.  Figure  12  of  Reference  22  shows  a near  conical 
cross  flow  separation  originating  at  approximately  16%  of  the  cones 
length  from  the  nose.  As  can  be  observed  from  the  oil  flow  traces  prior 
to  this  point,  the  flow  originating  at  or  near  the  nose  has  been  swept 
around  the  cone  to  the  lee  side  by  the  time  the  16%  point  has  been 
reached.  Since  the  flow  is  primarily  axial  on  the  lee  side  after  the 
cross  flow  separation  begins,  the  results  of  the  effective  blunting 
and  the  non-conical  region  near  the  nose  remain  at  or  near  the  lee 
centerline  as  the  flow  continues  along  the  cone. 

To  demonstrate  the  pressure  distribution  obtained  through  the 
shock  layer.  Figure  15  shows  static  pressure  in  the  6 direction  at  60° 
intervals  around  the  cone  for  the  a = 12°  case.  The  elevation  above 
the  cone  surface  is  given  in  radians  in  this  and  subsequent  computer 
generated  plots.  As  noted  in  Appendix  B,  the  damping  was  not  tailored 
in  the  0 direction  for  these  runs,  so  excess  smearing  of  the  shock  is 
evident  on  the  free  stream  side.  No  appreciable  harm  is  seen  to  result 
from  this  smearing. 

The  velocity  vector  plots  in  Figures  16  through  18  provide  an 
indication  of  cross-flow  streamline  patterns  and  illustrate  the  relative 
magnitude  of  cross-flow  velocity  for  all  but  the  lowest  momentum  regions. 
The  reverse  flow  region  is  thin  in  the  a = 12°  case  (Figure  16)  with 
the  vortical  singularity  (defined  for  this  study  as  the  cross  flow 
stagnation  point  toward  which  streamlines  converge)  occurring  near  the 
edge  of  the  viscous  layer  3.2°  above  the  cone  surface.  The  cross-flow 
separation  point  is  at  ())  = 163°  with  the  experimental  point  occurring  at 
(p  = 156°.  This  difference  is  attributable  to  lack  of  resolution  of  the 
reverse  flow  region  at  this  angle  of  attack  in  the  present  study,  as 
both  Lin  and  Rubin,  and  Lubard  and  Helliwell  obtained  a separation 
point  nearer  that  of  the  experimental  data.  This  was  supported  in  the 
present  study  by  the  fact  that  the  separation  points  in  the  higher  cx 
runs  (with  more  mesh  points  in  the  reverse  flow  regions)  were  virtually 
identical  to  those  of  experiment. 
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At  a = 20°  (Figure  17)  the  extent  of  the  reverse  flow  region  and 
the  position  of  the  cross-flow  separation  point  have  reached  limiting 
values,  as  no  substantial  changes  in  these  can  be  seen  in  the  a = 24° 
run  (Figure  18).  This  was  noted  for  the  separation  point  in  both  experi- 
ment and  theory  by  Lin  and  Rubin  (Reference  15).  However,  the  extent 
of  the  viscous  region  and  the  location  of  the  vortical  singularity  (as 
determined  in  the  present  study)  continue  to  change  as  a increases. 

The  vortical  singularity  is  near  the  edge  of  the  viscous  layer  for  all 
cases  in  the  present  study  in  which  it  is  lifted  off  the  cone  surface. 

Although  agreement  with  experiment  so  far  has  been  shown  to  be 
good,  perhaps  the  best  evaluation  of  the  validity  of  the  conical  flow 
assumption  for  engineering  solutions  can  be  made  through  a pictorial 
representation  of  the  flow  field  features.  Figures  19  and  20  compare 
the  overall  results  of  the  present  study  with  the  calculation  of  Lubard 
ad  Helliwell  and  Tracy's  data.  As  shown  in  Figure  19,  the  position  of 
the  bow  shock  and  viscous  layer  edge  of  the  present  study  for  this 
case  are  essentially  those  obtained  by  Lubard  and  Helliwell.  Both 
numerical  studies  result  in.  a shock  position  too  near  the  cone  and  a 
thinner  viscous  layer  at  this  a when  compared  with  experiment.  This 
is  considered  to  be  primarily  due  to  the  failure  to  account  for  nose 
tip  effects  as  noted  previously.  A weak  supersonic  region  is  present 
in  the  cross-flow  plane  and  is  shown  by  the  solid  lines  between  viscous 
layer  and  bow  shocK.  The' vortical  singularity  position  is  shown  by 
the  open  oval  symbol  on  the  lee  centerline.  The  structure  of  the 
a = 24°  case  is  more  interesting  (Figure  20).  Supersonic  cross-flow  is 
seen  to  exist  over  most  of  the  field  and  is  terminated  by  a complex 
sonic  line/shock  wave  The  position  of  this  line  is  nearer  the  lee 
centerline  than  shown  by  Tracy.  The  small  amplitude  oscillations 
present  on  the  imbedded  sonic  line/shock  as  it  nears  the  bow  shock  can 
be  seen  (through  careful  examination  of  Figure  18)  to  be  caused  by 
oscillations  (in  velocity)  propagating  through  the  supersonic  cross- 
flow  region  from  the  vicinity  of  the  bow  shock.  The  (J)  grid  spacing 
used  in  the  present  study  did  not  give  sufficient  resolution  to  deter- 
mine the  existence  of  the  lambda  portion  of  the  imbedded  shock,  as 
shown  by  Tracy.  However,  its  signature  may  be  surmised  in  the  hump 
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present  in  the  calculated  viscous  layer  edge  at  an  angle  of  approximately 
30°  from  the  lee  centerline.  The  curious  flat  top  present  in  the  cal- 
culated viscous  layer  edge  as  compared  to  experiment  is  again  attributed 
to  the  failure  of  the  numerics  to  model  non-conical  nose  effects. 

In  summary,  the  present  technique  is  seen  to  model  all  features 
present  in  the  experimental  study,  with  the  exception  of  the  non-conical 
nose  effects.  The  computer  times  are  not  excessive  since  they  are  based 
on  a very  stringent  convergence  criteria  and  represent  an  increase  in 
time  per  mesh  point  calculated  of  only  40%  over  a comparable  inviscid 
computer  program. 

2.  STETSON'S  DATA 

Stetson  (Reference  22)  conducted  af'  'experimental  study  using  a 5.6° 
half  angle  sharp  (0.001  in  nose  bluntness)  cone  in  H = 14.2  flow. 

The  wind  tunnel  working  fluid  was  air  with  total  pressure  of  160C  psia 
and  a total  temperature  of  2050°R.  The  resulting  free  stream  Re  is  given 
as  0.62  X 10^/ft.  Data  were  taken  at  angles  of  attack  of  0°  to  14°  in  2° 
increments  (plus  a run  at  5°).  Measurements  of  surface  pressure  were  made 
at  a number  of  stations  along  the  cone,  and  pitot  surveys  were  made 
perpendicular  to  the  free  stream  velocity  at  x/L  = 0.87  for  a = 10°. 

A calculation  was  made  at  the  following  conditions: 

M = 14.2  Re^  = 7.9  x 10^  (x/L  = 0.75, 

= 1600.0  psia  L = 15.37  in.) 

00 

T.  = 2050°R  T /T.  = 0.29 

t w t 

OO  00 

Pr  = 0.72 

The  finite  difference  mesh  contained  48  x 50  po^iits  with  the  circum- 
ferential step  size  being  4°  as  for  the  previous  calculations.  The 
computer  run  time  (with  convergence  criteria  as  before)  was  2.2  hours 
of  CDC  6600  time  for  2700  time  steps.  In  order  to  provide  comparison 
of  both  flow  field  and  surface  quantities,  plots  of  calculated  results 
versus  experimental  data  are  shown  for  surface  pressure  and  selected 
cuts  through  the  pitot  pressure  surveys.  Also,  a velocity  vector  plot 
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of  the  calculated  results  is  shown.  A plot  of  surface  pressure 
compared  with  Stetson's  data  comprises  Figure  21.  The  agreement  with 
experiment  is  quite  reasonable  with  the  maximum  difference  of  approx- 
imately 16%  (as  compared  with  local  static  pressure)  occurring  in  the 
vicinity  of  the  pressure  minimum.  Pressure  scale  expansion  in  Figure  21 
and  the  low  absolute  magnitude  of  pressure  in  this  region  make  the 
agreement  appear  less  adequate  than  actually  exists.  The  complete 
circumferential  pressure  distribution  is  qualitatively  comparable  to 
that  of  Figure  9.  The  surface  pressure  scale  of  Figure  21  is  equivalent 
to  0.0  to  0.6  on  Figure  6.  The  calculated  cross-flow  separation  point 
(26°  from  lee  centerline)  agrees  quite  well  with  the  experimental  value 
(27°  from  lee  centerline). 

Pitot  pressure  comparisons  at  three  stations  on  the  lee  side 
of  the  cone  are  given  in  Figures  22  through  24.  The  first  comparison  at 
(P  = 180°  (Figure  22)  shows  a large  difference  in  position  (approximately 
27%  at  the  bow  shock)  between  given  levels  of  pitot  pressure  for  the 
calculation  and  experiment.  Although  the  same  phenomenon  was  noted  in 
Tracy's  a = 24°  case  (Figure  14),  the  difference  was  small  at  the  edge 
of  the  viscous  layer  and  the  calculated  bow  shock  was  displaced  very 
little  from  the  experimental  position.  However,  in  the  present  case, 
the  bow  shock  is  displaced  nearly  as  much  as  the  edge  of  the  viscous 
layer.  This  displacement  can  be  seen  to  persist  (although  much  smaller 
in  magnitude)  through  = 160°  (Figure  23)  until  by  0 = 140°  (Figure  24), 
the  agreement  is  nearer  that  expected  from  previous  results.  In  all 
three  plots  of  pitot  pressure,  the  trends  as  shown  by  the  calculation 
are  essentially  correct.  The  remaining  question  then  concerns  the 
possible  reason  for  the  discrepancy  between  the  calculated  and  experi- 
mental viscous  layer  thickness. 

Several  possibilities  for  the  cause  of  the  viscous  layer  discrepancy 
should  be  examined.  It  should  first  be  noted  that  the  perfect  gas 
assumption  results  in  a somewhat  higher  Reynolds  number  than  is  quoted 
in  Reference  22  (0.83  x 10^  vs.  0.62  x 10^).  Since  this  would  result 
in  a somewhat  thinner  viscous  layer  for  the  calculation,  the  Reynolds 
number  was  lowered  20  and  the  run  was  reconverqed.  The  change  in 
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shock  position  on  the  leeward  side  was  a inaximuni  of  0.25°  (as  miqht 
have  been  expected  for  this  small  change  in  Re).  The  resulting  small 
change  in  shock  position  indicated  that  a change  in  Re  of  at  least  an 
order  of  magnitude  would  be  necessary  for  better  agreement.  A more 
likely  candidate  for  at  least  part  of  the  discrepancy  in  viscous  thick- 
ness is  the  effective  viscous  nose  blunting  present  in  the  experiment. 

This  is  clearly  shown  in  Figure  25  (taken  from  Reference  22)  and  would 
account  for  at  least  one-fourth  of  the  discrepancy  (assuming  no  change 
in  bow  shock  angle,  which  of  course  could  not  be  guaranteed). 

The  value  of  x foi"  this  calculation  was  3.22.  Without  further 
information  concerning  the  source  of  the  viscous  layer  discrepancy, 
this  higher  value  of  x must  be  viewed  as  an  indication  that  the  dis- 
crepancy will  be  encountered.  (According  to  Figure  6,  divergence  of 
the  present  technique  from  weak  interaction  theory  was  noted  at  x = 1-55 
for  Tracy's  conditions). 

To  illustrate  the  features  of  the  flow,  a plot  of  cross-flow 
velocity  vectors  comprises  Figure  26.  This  plot  illustrates  that  the 
lee  side  vortices  and  the  vortical  singularity  which  existed  in  Tracy's 
cases  exist  here  also.  The  vortical  singularity  occurs  near  the  edge 
of  the  viscous  layer  as  before. 

In  summary,  the  calculation  is  seen  to  model  the  features  and 
trends  present  in  Stetson's  experiniental  study  with  the  only  apparent 
difficulty  being  the  failure  to  correctly  predict  lee  side  viscous 
layer  thickness  and  shock  position.  The  maximum  difference  of  approx- 
imately 27*  occurred  at  the  lee  centerline  ((J  = 180°). 

3 i^cELDERRY'S  DATA 

McElderry  (Reference  31)  analyzed  experimental  data  taken  for  a 6' 
half  angle  sharp  cone  38.06  inches  in  length  in  M = 6.05  flow.  It  should 
be  noted  that  these  data,  also  appear  in  Reference  32  by  Rhudy  and  Baker. 

The  wind  tunnel  working  fluid  was  air  with  total  pressure  of  280.0  psia 
and  a total  temperature  of  850°R.  The  resulting  free  stream  Re  vjas 
given  as  5.0  x 10^/ft.  Data  were  taken  at  angles  of  attack  of  0°  to  12" 
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in  3°  increments.  Surface  pressure  measurements  were  made  at  several 
locations  along  the  cone  with  pitot  surveys  made  perpendicular  to  the 
cone  centerline  at  x/L  = 0.75  and  0.97.  Heat  transfer  data  were 
taken  on  a thin  walled  model  and  were  used  to  locate  the  boundary  layer- 
transition  region  on  the  cone.  References  to  transition  contained  in 
the  remainder  of  the  discussion  of  these  data  are  based  on  the  transition 
location  analysis  presented  in  Reference  31.  At  the  experimental  free 
stream  Reynolds  number  of  5 x 10^/ft  , transition  was  essentially  com- 
plete on  the  lee  side  by  x/L  - 0.3.  Since  the  surface  pressure  measure- 
ment nearest  the  tip  was  at  x/L  = 0.31,  nearly  all  lee  side  pressure 
data  were  taken  in  the  turbulent  flow  region. 


The  first  attempts  made  to  calculate  this  case  were  at  a Reynolds 
number  of  15.0  x 10^  which  co':  responds  to  the  fully  turbulent  flow 
region  near  the  base  of  tie  cone.  As  expected,  the  very  thin  boundary 
layer  present  at  this  Reynolds  number  created  extreme  difficulty  in  the 
vicinity  of  the  interaction  with  the  lee  side  imbedded  shock  wave. 

(The  flow  at  this  Reynolds  number,  when  considered  without  provision 
for  turbulence  modeling,  could  be  characterized  as  nearly  inviscid 
and  therefore  contains  some  of  the  difficulties  encountered  with  fhe 
inviscid  solution).  No  combination  of  damping  and  step  size  was  found 
which  would  lead  to  convergence  at  the  above  Reynolds  number.  The 
Reynolds  number  was  then  reduced  by  a factor  of  ten  to  1.5  x 10^  and 
the  run  was  successfully  comp''x.ted.  This  reduced  Reynolds  number 
corresponds  to  a .tation  on  the  cone  ju<'t  prior  to  the  onset  of  transi- 
tion (as  determined  in  the  experimental  study)  at  the  noted  angle  of 
attack  and  therefore  wa  in  a region  where  the  laminar  flow  equations 
were  still  valid. 

The  conditions  for  the  successful  run  were; 


M = 6.05 


= 280.0  psia 
Tj.  = 850  R 


Re  -1.5x10°  (x  4 in.  ) 

T /T.  = 0.61 

w t_ 

Pr  = 0.72 


32 


AFFDL-TR-76-139 


The  finite  difference  mesh  contained  48  x 50  points  with  the  c step 
size  being  4°  as  before.  Since  considerable  experimentation  with  time 
step  size  and  damping  factor  was  necessary  during  this  run,  the  computer 
time  of  2.3  hours  was  more  than  would  be  required  for  a rerun  of  this  or 
a similar  case. 

The  initial  behavior  of  this  solution  was  essentially  the  same  as 
that  for  the  higher  Reynolds  number  mentioned  above.  That  is,  the  lee 
side  imbedded  shock  is  strong  until  very  near  the  cone  surface  (due  to 
the  thin  viscous  layer).  This  resulted  in  a strong  compression  in  the 
vicinity  of  the  cone  surface  which  was  then  followed  by  an  overexpansion 
and  then  a recompression  to  the  cone  lee  centerline  point  (<p  - 180'  ). 
Although  this  appeared  to  be  a per'sistent  solution  to  tiie  governing 
equations  as  numerically  approximated,  the  overexpansion  became  strong 
enough  to  drop  pressure  (and  thereby  temperature)  below  zero  which  then 
destroyed  the  calculation.  Since  normal  stress  damping  was  found  to  be 
effective  in  damping  starting  transients  in  the  lee  side  region  for 
Stetson's  data,  it  was  also  tried  here.  After  some  experimentation,  it 
was  found  that  the  overexpansion  could  be  controlled  (and  eventually 
made  to  disappear  entirely)  by  a damping  factor  of  225.0  in  the  G 
direction  and  100.0  in  the  cj)  direction.  These  factors  are  scaled  by  the 
Reynolds  number  and  therefore  correspond  to  factors  of  63.0  and  28.0 
respectively  at  Tracy's  Reynolds  number.  In  contrast  to  the  use  of 
damping  to  control  starting  transients,  the  damping  in  this  case  could 
not  be  removed  as  the  overexpansion  would  then  reappear.  The  above 
damping  factors  were  then  maintained  until  the  solution  converged. 

Since  the  maintenance  of  normal  stress  damping  in  the  viscous  layer 
may  shift  the  density  distribution  through  the  layer  somewhat,  changes 
in  the  agreement  with  experiment  might  be  expected. 

Unfortunately,  experimental  measurements  were  not  made  far  enough 
forward  on  the  cone  to  provide  comparison  with  the  calculated  solution 
at  the  same  Reynolds  number.  The  first  surface  pressure  measurements 
were  at  x/L  = 0.31  at  which  point  transition  was  essentially  complete 
on  the  lee  side  at  12"  angle  of  attack.  The  dual  lee  side  separations 
associated  with  turbulent  flow  were  present  in  the  experiment  at  the 
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x/L  = 0.31  point  indicating  that  the  comparison  could  not  possibly  be 
exact  as  the  calculation  only  had  the  single  cross-flow  separation  which 
is  generally  associated  with  laminar  flow.  However,  this  calculation 
provides  a unigue  opportunity  to  examine  the  adeguacy  of  the  conical 
symmetry  assumption  for  turbulent  flow.  It  can  be  postulated  that  if 
the  conical  symmetry  assumption  is  to  be  valid  in  the  turbulent  regime, 
an  examination  of  the  flow  field  just  prior  to  the  start  of  transition 
would  differ  only  in  the  details  of  the  viscous  layer  when  compared 
conically  with  a measurement  far  downstream  of  transition.  To  test  this 
hypothesis,  a comparison  of  the  calculation  (x/L  = 0.10,  Re  = 1.5  x 10^) 
to  the  experimental  surface  pressure  (x/L  = 0.955,  Re  = 15.0  x 10^, 
Reference  31)  and  selected  pitot  surveys  (x/L  = 0.75,  Re  = 12.0  x 10^, 
Reference  32)  is  shown  in  Figures  27  through  29.  In  this  instance,  a 
conical  comparison  implies  examining  flow  quantities  at  two  distinct 
spherical  radii  (points  A and  B of  Figure  1)  on  each  of  a series  of  rays 
passing  through  the  vertex  of  the  cone.  If  the  flow  quantities  have 
nearly  the  same  vector  and  scalar  values  at  points  A and  B on  a majority 
of  the  rays,  then  the  flow  is  essentially  conical  in  character.  In  the 
present  comparison,  point  A corresponds  to  the  calculated  values  at  a 
point  on  the  cone  upstream  of  transition  and  point  B corresponds  to  the 
experimental  measurements  downstream  of  transition. 

The  experimental  surface  pressure  at  x/L  = 0.955  is  shown  conically 
compared  with  the  calculated  values  in  Figure  27.  The  agreement  is  quite 
good  in  both  trend  and  absolute  value  until  the  vicinity  of  the  experi- 
mental primary  separation  at  (!)  = 126°  is  reached.  From  this  point  on 
toward  the  lee  centerline,  the  absolute  value  of  the  pressure  shows 
varying  agreement  with  experiment.  The  calculated  separation  point 
(it  = 163°)  is  near  the  experimental  secondary  separation  point  (q  = 166"). 
To  place  the  conical  comparison  in  perspective,  the  agreement  of  the 
present  calculation  with  experiment  is  far  better  than  is  obtained  by 
the  inviscid  techniques  shown  in  Reference  31. 

Pitot  pressure  surveys  for  both  calculation  and  experiment  normal 
to  the  cone  centerline  are  presented  in  Fiqures  28  through  30.  The  first 
survey  is  at  ij)  = 180°  (Figure  28).  Both  experiment  and  calculation  are 
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characterized  by  very  high  speed  flow  prevailing  nearly  to  the  cone 
surface.  For  example,  the  first  two  mesh  points  above  the  surface  in 
the  calculation  have  local  Mach  numbers  of  5.0  and  5.9,  respectively. 

The  oscillations  occurring  between  1.0  and  2.0  inches  in  the  calculation 
are  numerical.  However,  the  apparent  oscillation  in  the  calculation 
occurring  between  0.0  and  1.0  inch  is  in  reality  the  solutions  reflection 
of  the  fact  that  a high  speed  "jet"  of  flow  occurs  (as  evidenced  by  the 
high  local  Mach  numbers  noted  above)  on  the  lee  centerline  under  the 
contra-rotating  vortices.  The  existence  of  this  jet  is  confirmed  by 
examination  of  Figure  31  which  shows  the  calculated  Mach  number  distri- 
bution along  lines  of  constant  $.  The  vortex  structure  is  seen  to  be 
centered  at  (p  170°.  The  local  Mach  number  clearly  increases  at  the 
first  two  mesh  points  above  the  surface  as  the  lee  centerline  ((p  = 180°) 
is  approached.  The  apparent  oscillation  on  the  0 = 180°  line  between 
0.0  and  1.0  inch  is  in  reality  an  indication  of  vortex  position.  This 
region  of  high  speed  flow  can  also  be  readily  seen  in  the  pitot  maps  in 
Reference  32  at  x/L  = 0.97.  However,  the  resolution  of  the  lines  of 
constant  pitot  pressure  is  insufficient  for  the  x/L  = 0.75  surveys  in 
Reference  32  to  be  able  to  plot  with  certainty  this  "jet"  in  the  experi- 
mental data  on  the  graph  shown  here.  The  pitot  survey  at  (f  = 160° 

(Figure  29)  shows  much  the  same  trend  as  the  ({>  = 180°  survey  except  for 
lower  momentum  flow  near  the  surface  and  a slightly  thinner  viscous 
layer  in  the  calculation  than  in  the  experiment.  At  (p  = 140°  (Figure  30) 
the  conical  agreement  of  the  calculation  with  experiment  is  very  good. 

The  global  conical  agreement  with  experiment  can  be  examined  through 
the  characteristics  of  the  cross-flow  velocity  vector  plot  in  Figure  32. 

One  immediately  obvious  feature  of  this  plot  is  that  the  vertical  singu- 
larity is  not  lifted  off  the  surface  even  though  the  angle  of  attack  is 
twice  the  cone  half  angle.  This  is  in  reality  evidence  of  the  fact  that 
the  viscous  layer  instead  of  having  a large  hump  on  the  lee  side  now  has 
a "cusp"  at  the  lee  centerline  with  a nearly  sharp  local  minimum  in 
viscous  layer  thickness.  The  experimentalist  have  observed  this  behavior 
of  the  viscous  layer  edge  and  it  has  been  noted  in  vapor  screens  and 
pitot  surveys  taken  in  turbulent  flow  on  cones  in  References  31,  33,  and  34 
(see  Figure  33  taken  from  Reference  33).  Also,  comparing  the  extent  of 
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the  lee  side  vortices  present  in  the  laminar  reqion  calculation  with  the 
pitot  maps  in  References  31  and  32,  it  is  apparent  that  the  essential 
size  and  character  of  the  primary  vortices  are  maintained  through 
transition  and  beyond. 

The  differences  in  the  characteristics  of  the  viscous  layer  and 
surface  pressure  distribution  with  increased  Reynolds  numbers  are 
graphically  and  picti^-ially  shown  in  Figure  33  which  is  taken  from 
Reference  33.  The  "cusp"  viscous  layer  structure  shown  in  Figure  33 
for  Rainbird's  experimental  results  has  been  considered  to  be  primarily  a 
turbulent  flow  regime  phenomenon  on  conical  bodies.  However,  the  calcu- 
lation made  at  McElderry's  conditions  in  the  present  study  demonstrates 
that  the  "cusp"  is  present  in  the  viscous  layer  upstream  of  transition. 
The  viscous  layer  structure  of  Rainbird  shown  in  Figure  33  is  therefore 
a phenomenon  occurring  in  the  laminar  flow  regime  which  persists  into 
the  turbulent  regime. 

In  summary,  the  solution  at  a location  on  the  cone  just  upstream  of 
the  boundary  layer  transition  location  is  shown  to  exhibit  features 
previously  associated  with  turbulent  flows  over  conical  bodies.  The 
agreement  of  this  solution  with  experiment  is  amazingly  good  when  it 
is  compared  conically  with  experimental  measurements  made  far  down- 
stream in  the  turbulent  flow  regime.  The  adequacy  of  the  conical 
symmetry  assumption  for  turbulent  flow  is  therefore  assured  if  a 
reasonable  means  can  be  found  to  model  viscous  layer  turbulence. 
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SECTION  V 


CONCLUSIONS 


Inviscid  computations  for  hypersonic  cone  flow  have 

unstable  at  angles  of  attack  greater  than  twice  the  cone 

In  addition,  previously  obtained  solutions  including  the 

viscosity  have  required  input  of  experimental ly  obtained 

have  not  been  demonstrated  for  a/R  >2.0. 

c 


proven  to  be 
half  angle, 
effects  of 
data  and/or 


In  order  to  resolve  these  difficulties,  solutions  have  been 
demonstrated  for  the  complete  flow  field  around  sharp  cones  including 
viscous  effects.  These  solutions  were  obtained  by  applying  a conical 
symmetry  assumption  to  the  complete  Navier-Stokes  equations  and  numerically 
integrating  the  resulting  equation  set  by  use  of  MacCormack's  finite 
difference  scheme.  Integrations  were  performed  at  selected  sets  of 
boundary  conditions  for  which  previously  published  experimental  data 
were  available.  The  solutions  obtained  were  for  a Mach  number  range  of 

r r 

6.05  to  14.2,  a local  Reynolds  number  range  of  0.4  x 10°  to  1.5  x 10  , 
and  angles  of  attack  from  a/9^  = 0.0  to  2.4.  Stability  of  the  solutions 
was  demonstrated  through  use  of  a stringent  convergence  requirement. 

A technique  (normal  stress  damping)  was  developed  and  used  to  reduce  or 
eliminate  the  numerical  oscillations  which  developed  during  the  inte- 
gration in  the  vicinity  of  strong  shock  waves.  The  results  of  the 
integration  were  compared  with  the  experimental  studies  of  Tracy, 

Stetson,  and  McElderry.  The  Following  points  are  presented  concerning 
the  conduct  of  the  integration  and  the  comparison  with  experiment: 


(1)  Agreement  with  the  experimental  data  of  Tracy  was  generally 
good.  The  major  features  of  the  experiment  (shock  waves,  sonic  lines,  and 
viscous  layers)  were  present  and  essentially  correct  in  location  and 
magnitude  in  the  results  of  the  integration.  The  primary  discrepancies 
noted  in  the  comparison  with  Tracy's  data  were  an  approximate  5“''  dif- 
ference in  surface  pressure  near  the  windward  centerline  and  an  error  in 
pitot  pressure  in  a small  region  near  the  lee  centerline  for  the 
i/0^  2.0  cases.  The  latter  discrepancy  was  concluded  to  be  due  to 
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failure  of  the  present  technique  to  account  for  non-conical  nose  effects. 
Agreement  of  the  results  of  the  integration  with  the  experimental  data  of 
Stetson  was  less  adequate  than  was  shown  to  exist  with  the  data  of  Tracy. 

The  agreement  with  Stetson's  surface  pressure  measurements  was  quite 
reasonable.  However,  the  position  of  the  lee  side  shock  wave  in  the 
results  was  27%  nearer  the  cone  surface  at  the  lee  centerline  than  in 
the  experiment.  This  discrepancy  was  concluded  to  be  primarily  caused 
by  inability  of  the  present  technique  to  model  the  effective  nose 
blunting  existing  in  the  experiment  which  resulted  from  viscous  effects 
at  this  high  Mach  number  (14.2).  The  apparent  effect  of  this  nose 
blunting  was  much  greater  in  extent  in  Stetson's  case  (hypersonic  similarity 
parameter,  x - 3.0)  than  in  Tracy's  cases  (x  ~ 1.0). 

(2)  The  present  technique  removes  the  angle  of  attack  limitations 
present  in  inviscid  calculations.  The  instabilities  associated  with  the 
point  where  the  high  strength  lee  side  imbedded  shock  wave  reaches  the 
surface  are  removed  by  the  inclusion  of  viscous  effects  in  the  physical 
modeling  of  the  flow  field. 

(3)  The  adequacy  of  the  conical  symmetry  assumption  is  indicated 
for  the  turbulent  regime  on  conical  bodies.  A solution  obtained  just 
prior  to  the  beginning  of  boundary  layer  transition  at  the  experimental 
conditions  of  McElderry  agrees  well  with  experiment  when  conically 
projected  far  into  the  turbulent  regime. 

(4)  Normal  stress  damping  is  shown  to  provide  a physically  based 
means  to  provide  the  necessary  control  of  spurious  numerical  oscillations 
around  shock  waves  without  additional  computer  time.  It  is  also  very 
effective  for  control  of  overexpansion  and  for  control  of  starting 
transients  due  to  ill-suited  initial  conditions. 

In  summary,  the  present  technique  was  shown  to  be  a viable  means  of 
calculating  the  flow  over  cones  at  angle  of  attack  including  viscous 
effects.  Agreement  with  experiment  was  quite  good  for  x £ 1 with  some 
discrepancy  encountered  in  the  lee  side  viscous  layer  as  x increased. 
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Angle  of  attack  limitations  encountered  with  the  inviscid  equations 
were  removed.  This  technique  should  therefore  be  considered  as  an 
alternative  to  the  use  of  the  inviscid  flow  equations  for  future 
calculations  of  flow  over  conical  bodies. 
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SECTION  VI 
RECOMMENDATIONS 

Three  primary  directions  are  seen  for  extension  of  the  work  presented 
in  this  report.  The  first  is  the  confirmation  of  the  technique  for 
more  general  conical  configurations,  such  as  delta  wings.  The  second  is 
the  extension  of  the  Reynolds  number  range  for  which  the  calculations  can 
be  accomplished  through  inclusion  of  a turbulence  model  in  the  formulation 
of  the  equation  set.  The  third  is  the  investigation  of  the  adequacy  of 
the  conical  symmetry  assumption  for  lower  supersonic  Mach  numbers. 

A further  item  pointed  out  by  this  study  is  the  need  for  re- 
examination by  the  experimentalist  of  the  techniques  used  for  collection 
of  experimental  data.  As  the  numerical  computation  capability  improves, 
heretofore  undetected  discrepancies  in  data  collection  techniques  may  be 
pointed  out.  An  example  is  the  underestimation  of  surface  pressure  on 
the  windward  side  of  the  cone  in  the  present  study.  This  underestimation 
in  past  inviscid  studies  was  attributed  to  failure  to  include  viscous 
effects.  However,  two  studies  (the  present  study  and  Reference  16)  are 
available  which  exhibit  this  underestimation  even  though  viscous  effects 
were  included  in  the  calculated  results.  A third  study  (Reference  14) 
matched  the  windward  surface  pressure  for  Tracy's  data  only  at  a much 
lower  Reynolds  number  than  existed  in  the  data.  This  is  considered  to 
be  a clear  indication  that  the  details  of  the  use  of  surface  pressure 
taps  where  thin  hypersonic  boundary  layers  exist  should  be  further 
studied.  This  also  points  out  the  need  for  cooperation  and  collaboration 
between  the  experimentalist  and  the  numericist  in  the  effort  to 
improve  overall  aerodynamic  prediction  techniques. 
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Figure  6.  Comparison  of  Calculated  Cone  Surface  Pressure  with 

that  Obtained  from  Hypersonic  Weak  Interaction  Theory 
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Figure  7.  Comparison  of  Calculated  and  Experimental  Surface 
Pressure  at  Tracy's  Conditions,  a = 8° 
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Figure  8.  Comparison  of  Calculated  and  Experimental  Surface 
Pressure  at  Tracy's  Conditions,  a = 12° 
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Figure  9.  Comparison  of  Calculated  and  Experimental  Surface 
Pressure  at  Tracy's  Conditions,  a = 20° 
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Figure  10.  Comparison  of  Calculated  and  Experimental  Surface 
Pressure  at  Tracy's  Conditions,  a = 24° 
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Figure  16.  Cross-Flow  Velocity  Vectors  Plotted  as  Seen  on  Spherical 
Surface  at  Tracy's  Conditions,  a = 12° 
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! Figure  18.  Cross-Flow  Velocity  Vectors  Plotted  as  Seen  on  Spherical  I 

I Surface  at  Tracy's  Conditions,  a = 24°  j 
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Figure  19.  Comparison  of  Calculated  and  Experimental  Flow  Field 
Features  at  Tracy's  Conditions,  a = 12" 
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Figure  21.  Comparison  of  Calculated  and  E 
at  Stetson's  Conditions 
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Figure  22.  Comparison  of  Pitot  Pressure  Survey  at  Constant  (}i  at 
Stetson's  Conditions,  <})  = 180° 
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Figure  23.  Comparison  of  Pitot  Pressure  Survey  at  Constant  (p  at 
Stetson's  Conditions,  A = 160° 
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1 Figure  24.  Comparison  of  Pitot.  Pressure  Survey  at  Constant 

j Stetson's  Conditions,  o = 140° 
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Figure  28.  Conical  Comparison  of  Calculated  (Re  = 1.5  x 10^)  and 

6 ^ 

Experimental  (Re  = 15.0  x 10  ) Pitot  Pressure  Survey 
at  Constant  4)  at  McElderry's  Conditions,  4 = 180° 
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Figure  29.  Conical  Comparison  of  Calculated  (Re  = 1.5  x 10^)  and 

6 ^ 

Experimental  (Re^  = 15.0  x 10°)  Pitot  Pressure  Survey 
at  Constant  4)  at  McElderry's  Conditions,  (|)  = 160° 


AFFDL-TR-76-139 


■ / 

✓ 

* 

^ Jt 


/ / 

/ / ^ 

/ ^ ^ 
/ / r 


iffss;5s:^ 

ili^ 


0:55:^ 


ELEVATION  IRQ.  ) 


Figure  32.  Cross-Flow  Velocity  Vectors  Plotted  as  Seen  on  Spherical 
Surface  at  McElderry's  Conditions 
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All  Heasurements  Have  Been  Normalized  by  the  Corresponding 
Zero  Angle-o£-Attack  Value. 


Figure  33.  Comparison  of  Experimental  Laminar  and  Turbulent  Viscous 
Layer  Features  (Taken  from  Reference  33) 


76 


P/?1N 


AFFDL-TR-76-139 


3 


Figure  34.  Demonstration  of  the  Effect  of  Normal  Stress  Damping  on 
Calculated  Static  Pressure  Along  6,  a = 0° 
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APPENDIX  A 
BOUNDARY  CONDITIONS 

The  solution  of  the  Navier-Stokes  equations  by  numerical  methods  is 
in  reality  a boundary  value  problem  subject  to  the  same  requirements  and 
constraints  as  the  solution  of  the  set  by  analytic  techniques.  However, 
the  mathematical  theory  concerning  existence  and  uniqueness  of  solutions 
for  this  equation  set  has  not  been  resolved  (Reference  35).  This  statement 
is  especially  true  for  compressible  flows.  The  only  course  then  open  to 
the  engineer  is  to  presume  a well  posed  problem  exists  when  the  known 
physically  derived  auxiliary  conditions  are  imposed.  If  the  integration  of 
the  equation  set  subject  to  these  conditions  yields  a steady  state  solution 
(which  of  course  is  not  guaranteed),  then  experimental  results  for  these 
same  conditions  must  be  used  to  verify  that  the  solution  obtained  is 
indeed  the  correct  one.  When  this  is  done  for  several  related  sets  of 
conditions  with  good  results,  confidence  in  the  adequacy  of  the  solution 
technique  for  use  where  experimental  results  are  not  available  is  increased. 
However,  the  user  of  the  technique  for  this  purpose  must  always  take  care 
to  examine  the  solutions  obtained  for  nonphysical  features,  as  it  is  never 
possible  to  completely  define  initially  the  limits  of  usefulness  of  an 
equation  set  in  any  meaningful  sense.  The  above  procedure  was  followed  in 
this  study.  The  rest  of  this  appendix  will  present  the  auxiliary  condi- 
tions used  and  the  manner  in  which  they  were  imposed. 

Examination  of  Figure  2 reveals  that  the  finite  difference  mesh 
is  a spherical  surface  with  the  inner  edge  bounded  by  the  cone  and 
the  outer  edge  free.  The  mesh  shown  is  360°  circumferentially  in 
extent.  Since  two  knowr  lines  of  symmetry  exist  at  (J)  = 0°  and  (|)  = 180°, 
the  extent  of  the  surface  can  be  readily  reduced  by  one-half.  This 
reduction  can  be  made  without  increase  in  the  number  of  boundary  values 
necessary  to  the  calculation  by  specifying  an  additional  line  of  mesh 
points  in  the  9 direction  and  at  the  required  value  of  <J).  The  results 
of  each  sweep  during  the  integration  are  then  reflected  symmetrically 
to  these  lines  about  the  lines  of  symmetry.  This  causes  the  inte- 
gration (at  the  line  of  symmetry)  to  perceive  that  a mirror  image  of 
the  solution  is  evolving  on  the  other  180°  of  the  mesh  (as  it  would 
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be  if  the  sweep  were  actually  360°).  This  procedure  therefore  re- 
quires no  a priori  knowledge  of  the  flow  quantities  at  the  lines  of 
symmetry  and  implies  that  the  following  conditions  are  being  imposed 
at  these  lines  of  symmetry: 


at  ({)  = 0° 

and  0 = 

180° 

M = 

30 

' Note  that  these  two  conditions 

0 1 

are  imposed  only  at  ^n_e  line 
of  symmetry.  They  are  then 

94)  ^ 

9p  _ 

0 1 

derivable  from  the  governing 

30 

[ equations  at  the  other. 

w = 0 

^ = 

0 

30 

Two  boundaries  therefore  remain  which  must  be  dealt  with. 


Examination  of  the  six  equations  (the  Navier-Stokes  equations,  the 
continuity  and  energy  equations,  plus  the  equation  of  state)  as  modified 
reveals  that  all  three  velocity  components  and  temperature  appear  as 
second  derivatives  in  0 and  4>  with  pressure  and  density  appearing  as 
first  derivatives  in  6 and  <}).  A solution  should  then  be  obtainable 
through  the  specification  of  two  boundary  conditions  for  each  of  the 
velocity  components  and  temperature  and  one  each  for  pressure  and  density. 
(In  reality,  the  algebraic  equation  of  state  eliminates  the  requirement 
for  specifying  both  pressure  and  density.)  This  was  done  in  the  present 
study  in  the  following  manner: 


at  the  cone  surface  0 = 0^, 


= 0°-vl80° 


u = V = w = 0 

T = T (<f) 
w 

at  the  outer  boundary  (physically,  the  free  stream) 
0 = 0 outer.  Ip  = 0°  >-180° 

u = u^  ((i)  p = p^ 

V = v^  (()))  p = 

VI  = VI  ( Ip)  T = T 


a uniform  free 
stream  is  assumed 
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Note  that  although  the  quantities  are  specified  along  a line  of  constant 
0,  the  Ip  distribution  is  completely  specified  also.  It  should  also  be 
noted  that  the  total  internal  energy  was  used  in  place  of  the  temperature 

U2 

as  a dependent  variable.  The  connection  is  simply  e = T + -2-  . 

These  conditions  are  now  sufficient  to  allow  solution  of  the  problem. 


Since  pressure  and  density  at  the  cone  surface  remain  free  (and 
indeed  must  do  so  to  prevent  overspecification  of  the  boundary  conditions), 
they  must  be  determined  during  the  integration  process.  This  was 
done  in  the  present  study  by  analytically  evaluating  the  0 momentum 
equation  ’t  the  cone  surface  conditions.  The  resulting  equation  is 
non-dimensional ized  form  solved  for  is: 


3P  - 

1 

! ^ 

r COS0 

+ sine  (. 

30 

k 

sin0 

1 3 

L Re 

8 

sin0 

^ + 

4 sin6 

3^v  2 3 

3 

Re 

36  ^ 

3 Re 

^ ■ 3 30 

( 


1 9w  \ 


'Re  3(J)' 


+ ^(1_  1.)  ! 
3(f)^Re  30^  I 


This  equation  is  approximated  by  one  sided  second-order  accurate  finite 
differences  about  the  surface  mesh  point  of  interest.  The  resulting 
difference  equation  is  solved  for  the  surface  pressure,  the  only  quantity 
in  the  equation  which  has  not  been  updated  when  the  equation  is  applied 
after  a sweep  through  the  mesh  of  either  the  predictor  or  the  corrector. 
This  equation  can  be  represented  by: 


.n 


- 2A0 
3k  sin0 


f 


(u,v,w) 


+ 


4 i 

3 Pj+l,n 


1 

3 


P 


i 

j+2,n 


with  j = 1,  and  i and  m indicating  time  level  and  $ position  in  the 
mesh  respectively.  The  density  at  the  surface  is  found  via  the  equation 
of  state  and  all  flow  quantities  are  then  known  or  have  been  determined 
at  the  cone  surface. 

To  complete  the  specification  of  the  auxiliary  conditions,  initial 
values  of  the  flow  quantities  are  required  at  all  points  in  the  mesh  in 
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order  to  begin  the  calculation.  In  the  present  study,  this  was  done  by 
setting  the  mesh  to  free  stream  conditions  including  the  mesh  point  above 
the  cone  surface  with  true  surface  boundary  conditions  plus  free  stream 
pressure  and  density  imposed  at  the  cone  surface.  The  calculation  was 
then  begun  which  resulted  in  a so  called  "impulsive"  start  up.  (i.e., 
the  cone  is  brought  from  rest  to  free  stream  velocity  instantaneously). 
The  outer  and  inner  boundary  conditions  were  maintained  throughout  the 
calculation. 

In  summary,  a means  of  reducing  the  field  of  calculation  by  half 
through  symmetry  has  been  discussed.  The  requirements  for  boundary 
conditions  were  shown  and  the  manner  in  which  these  requirements  were 
satisfied  in  the  present  study  was  presented.  The  use  of  the  0 momentum 
equation  to  obtain  the  cone  surface  pressure  was  demonstrated  and  the 
initial  conditions  utilized  were  also  discussed. 
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APPENDIX  B 

NORMAL  STRESS  DAMPING 

Previous  studies  (References  13,  25)  using  shock  capturing  finite 
difference  techniques  have  encountered  oscillations  in  the  vicinity  of 
strong  shock  vjaves  which  can  cause  solution  instabilities.  To  overcome 
this  difficulty,  artificial  terms  have  been  added  to  the  governing 
equations  to  provide  necessary  damping.  In  the  present  study,  using 
the  viscous  equations  and  shock  capturing,  oscillations  were  again 
encountered.  This  indicated  that  the  natural  viscous  terms  are  inade- 
quate to  overcome  the  series  truncation  error  at  the  existing  mesh 
spacing.  To  damp  these  oscillations  without  additional  programming 
complexity,  the  normal  stress  terms  (which  are  in  general  large  only  near 
shock  waves)  were  altered  by  increasing  the  second  coefficient  of 
viscosity  - 2/3u).  This  resulted  in  elimination  or  reduction  of  the 
oscillations  where  desired.  In  addition  to  improving  the  shock  structure, 
it  was  found  that  normal  stress  damping  was  extremely  effective  in 
removing  instabilities  caused  by  starting  transients  which  result  from 
ill-suited  initial  conditions.  The  details  of  implementation  and 
conditions  for  proper  use  of  this  technique  are  given  below. 

As  shown  in  Equation  4,  the  normal  stress  terms  contain: 

% = ^Re  = - iRe  [ ^ ^ If  ] 

Normal  stress  damping  results  when  2u  fin0  is  removed  and  the 
remainder  of  the  term  is  multiplied  by  an  empirically  determined  factor 
B,  where  P,  is  negative  (i.e.  8o^).  (A  positive  6 does  not  produce  the 
damping  effect  and  makes  the  calculation  unstable  except  near  B = 1). 

The  damping  terms  occur  both  in  the  energy  equation  and  in  the  H'  matrix 
resulting  from  the  coordinate  transformation.  It  was  found  that  the 
most  accurate  solutions  resulted  from  use  of  the  damping  multiplier  in 
the  derivative  term  only  (f^  or  d^)  of  the  momentum  equation  most 
nearly  normal  to  the  shock  wave  for  which  damping  is  desired.  Inclusion 
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of  the  multiplier  in  the  R'  matrix  or  the  energy  equation  produced  some 
displacement  of  the  solution.  Also,  it  should  be  noted  that  in  the  steady 
limit, 

= f(Vv)  = f(^^)  (17) 

(through  use  of  the  continuity  equation).  This  implies  that  normal 
stress  damping  cannot  be  used  in  the  vicinity  of  compressible  flow 
boundary  layers  without  regard  for  possible  changes  in  the  density 
gradients. 

An  example  of  the  use  of  normal  stress  damping  is  given  in  Figure 
34  for  axisymmetric  cone  flow  as  calculated  by  the  technique  of  this 
report.  The  plot  is  a comparison  of  the  calculated  static  pressure  in 
the  0 direction  at  Tracy's  free  stream  conditions  with  and  without  the 
inclusion  of  normal  stress  damping.  In  this  case  the  damping  is 
tailored  in  the  0 direction  by: 

Ps  P 

B = 1 - (^)  A1  - (^)  A2  (18) 


where  the  pressure  ratio  multipliers  were  empirically  determined  to  be 
A1  = 9 and  A2  = 2.  The  use  of  the  damping  multiplier  (Figure  34)  resulted 
in  the  reduction  of  the  numerical  oscillations  near  the  shock  wave. 
Although  a small  offset  in  pressure  occurs  outside  the  boundary  layer, 
the  surface  pressure  and  shock  position  do  not  change  appreciably 
when  damping  is  applied. 


For  Tracy's  cases,  the  damping  multiplier  B was  unity  over  the  lee 
40°  ((})  = 140°  - 180°)  of  the  flow,  where  viscous  effects  predominate. 

It  was  then  varied  smoothly  around  the  cone  from  <p  = 140°  in  accord 
with  the  equation:  , 


6 


*'sU  - ^sU  = 140° 

’’sU  = 0°  - ^s'cl'  = 140° 


(19) 


where  s indicates  body  surface  values  and  B'  was  empirically  varied  from 
15  (a  = 8°)  to  130  (a  = 24°).  B was  constant  in  the  0 direction. 
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For  Stetson's  case,  the  very  strong  lee  side  expansion  during  start 
up  made  it  necessary  to  use  normal  stress  damping  there.  However,  the 
damping  in  this  region  was  removed  after  the  solution  development  was 
sufficient  for  the  expansion  to  become  less  strong.  For  the  rest  of  the 
field  (<p  = 0°->140‘’),  a more  sophisticated  tailoring  of  the  damping 
multiplier  was  used: 


^ ^ = 1 - 
6 ,(() 

P L -P 

S l(J)  s 

1 

0 

O, 

II 

-e- 

1/2 

'"l.t  = 0°  1 

^14.  = 0 

° - ^s|<()  = 140° 

O 

O 

II 

Q. 

J 

'"U,e 

co2 

p 

00 

1 

value  of  col  was  1.6  and  of  co2  was  1.1 

This 

equation  was 

col 


(20) 


on  experiments  with  the  use  of  normal  stress  damping  which  showed  that 
the  effects  could  be  improved  by  tailoring  according  to  the  local  value 
of  static  pressure.  An  improvement  was  noted  in  the  damping  effect,  but 
it  was  found  that  this  equation  did  not  always  increase  the  magnitude 
of  6 quickly  enough  during  initial  start  up.  This  was  overcome  by  using 
a high  initial  magnitude  of  6 for  approximately  the  first  100  time  steps 
of  the  run.  The  above  equation  was  used  for  McElderry's  case  also. 
However,  the  value  of  B necessary  on  the  lee  side  was  more  negative  than 
the  output  of  the  above  equation.  The  high  initial  magnitude  was  there- 
fore used  everywhere  in  the  flow  field. 


In  order  that  others  who  use  normal  stress  damping  may  benefit  from 
experience  gained  in  the  present  study,  the  following  observations  and 
rules  concerning  this  technique  are  offered: 

(1)  The  most  valuable  and  recommended  use  of  normal  stress  damping 
is  to  reduce  starting  transients  from  ill-suited  initial  conditions.  No 
permanent  history  of  the  use  of  damping  for  this  purpose  is  seen  in  the 
flow  after  the  removal  of  the  damping  and  subsequent  convergence. 
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(2)  Normal  stress  damping  has  been  used  to  allow  capturing  of  shock 
waves  with  P2/P1  as  high  as  32.  The  only  practical  limit  to  shock  strength 
which  can  be  captured  through  use  of  this  technique  appears  to  be  number 

of  mesh  points  available  for  the  smearing  to  take  place  over. 

(3)  Tailoring  of  the  damping  based  on  local  static  pressure  was 
found  to  give  reasonable  results.  Tailoring  based  on  other  quantities 
(such  as  temperature)  gave  unnecessary  offset  in  the  viscous  layer. 

(4)  Large  shifts  of  surface  pressure  or  of  static  pressure  behind 
shocks  when  damping  is  applied  usually  indicate  that  the  magnitude  of  6 
is  too  high  in  that  locale. 

(5)  The  best  results  are  obtained  behind  the  shock  wave  when  some 
oscillations  are  allowed  to  remain  in  the  free  stream. 

(6)  Care  must  be  taken  to  ensure  that  the  viscous  stability  limit 
is  satisfied  with  g included  in  the  manner  she  " in  Appendix  C. 

(7)  The  items  noted  below  are  observe,  luse  undue  displacement 
of  the  solution  at  high  values  of  g: 

(a)  Inclusion  of  the  g multiplier  in  terms  resulting  from 
coordinate  transformations  (i.e.,  the  matrix). 

(b)  Inclusion  of  g in  the  energy  equation. 

The  previous  observations  and  rules  should  not  be  considered  as  axiomatic 
but  are  offered  as  guides  for  the  use  of  this  tp'-h-  oue  in  other  com- 
pressible flow  problems. 

Normal  stress  damping  is  a physically  based  technique  which  has 
proven  beneficial  in  damping  out  starting  transients  and  smoothing 
numerical  oscillations  in  the  vicinity  of  shock  waves.  The  B contouring 
equations  and  levels  used  for  the  various  cases  were  presented.  Rules 
and  observations  were  also  given  to  guide  future  use  of  this  technique. 
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APPENDIX  C 
STABILITY  CRITERIA 


In  order  to  maximize  computational  efficiency,  criteria  should  be 
used  to  determine  the  maximum  permissible  time  step  (At)  as  the 
integration  progresses.  As  noted  in  Section  III,  criteria  were  not  used 
during  the  calculations  presented  elsewhere  in  this  report.  Instead, 

At  was  fixed  during  each  complete  integration  at  Tracy's  conditions 
and  was  fixed  for  the  later  parts  of  the  integrations  at  Stetson's  and 
McElderry's  conditions.  The  value  at  which  At  was  fixed  was  at  first 
determined  by  increasing  At  until  instability  was  encountered  during 
the  starting  transient  phase  of  the  run.  The  run  could  then  be  continued 
to  convergence  by  reducing  At  10%  to  15%  from  the  unstable  value.  As 
experience  was  gained  with  the  technique  (and  more  was  learned  about  the 
effect  of  normal  stress  damping  on  stability),  it  was  found  that  for 
Tracy's  conditions,  stable  values  of  At  could  be  set  initially  without 
resorting  to  trial  and  error.  These  values  were  in  general  0.6A6  to 
O.7A0.  The  above  procedure  did  not,  however,  allow  for  any  possible 
increase  in  allowable  At  after  the  large  transient  phase  of  the  calcu- 
lation. Difficulty  in  determining  allowable  At  was  also  encountered 
for  Stetson's  and  McElderry's  conditions  due  to  use  of  normal  stress 
damping  in  predominately  viscous  regions.  To  correct  these  two 
deficiencies,  stability  criteria  for  use  during  the  integration  were 
postulated  and  confirmed  by  a repeat  of  previous  calculations.  The 
remainder  of  this  appendix  will  present  this  procedure  and  the  observed 
changes  in  computational  efficiency. 

A complete  stability  analysis  for  MacCormack's  method  as  applied 
to  the  Navier-Stokes  equations  has  not  been  accomplished.  However, 
useful  estimates  of  allowable  step  sizes  can  be  made  through  the  use  of 
linear  theory  and  physical  arguments.  Claiming  that  the  physical 
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propagation  of  information  (i.e.  fluid  signals)  cannot  outstrip  the 
numerical  propagation  results  in  the  following  Courant-Friedrich-Lewy  (CFL) 
condition  (Reference  36); 


At 


cfl  - 


A0 


sin9  A(Ji 


_1_,_0_ 

lAe*^  ___ 

' A6(sin0  Act) 


+ (sine  A(|)^] 


(21) 


This  expression  essentially  determines  allowable  At  in  regions  where 
viscous  effects  are  negligible.  The  remaining  task  is  to  determine  the 
effects  on  stability  of  the  stress  terms  and  of  normal  stress  damping. 


The  maximum  stable  At  in  viscous  regions  is  normally  a function  of 
the  smallest  physical  step  size  there.  In  the  present  study  A6«A<t 
for  all  cases  calculated.  This  implied  that  the  stress  terms  in  the  6 
coordinate  direction  were  the  main  determinant  of  allowable  At  and 
therefore  that  the  stability  analysis  due  to  the  stress  terms  was  essen- 
tially a one-dimensional  problem.  In  the  case  of  finite  difference 
operator  time  splitting,  the  stability  analysis  also  has  been  shown  to 
be  a one-dimensional  problem  (Reference  25).  An  analysis  similar  to 
that  of  Reference  25  was  therefore  accomplished  and  is  presented  here. 

Equation  22  is  first  written  in  nonconservation  form: 


[L][Z]  ^ + [M][Z]  Q + [N][Z]  +[W][Z]  + [Q][Z] 


❖4> 


- [s][zl  - [x][z]  .[H'l  = 0 


(22) 


where  [Z] 


P 

u 

V 

w 

p 


and  the  subscripts  indicate  differentiation. 


The  remaining  matrices  will  be  defined  as  they  are  used.  The  final 
equation  results  when  [L]~^  is  found  and  the  coefficient  matrices  [L]'^  [M], 
etc.  are  formed.  By  declaring  the  coefficient  matrices  to  be  locally 
constant,  it  can  be  shown  that  (References  37,  38)  when  each  dimension  is 
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considered  separately  1.0  (inviscid  effect)  and  ^ 1 ^ (viscous 

effect).  The  symbol  v is  in  each  case  the  maximum  eigenvalue  of  the 
appropriate  coefficient  matrix  (for  example  [L]'^  [M]  ).  Equation  21 
accounts  for  the  inviscid  effect  on  stability  so  that  it  is  only 
necessary  to  consider  the  second  inequality  and  obtain  values  of  v for 
the  coefficient  matrices  of  the  second  derivative  terms  involving  0 
in  Equation  22. 

The  matrix  [L]”^  [W]  is  as  follows: 


[L]'^  [W] 


0 

iA 

2 + y 
pRe 


2p2RePr 


pRe  Pr 


The  maximum  eigenvalues  of  this  matrix  are  y/(pRe  Pr)  and  (2  + ^)  /(pRe). 
The  second  eigenvalue  is  included  due  to  the  undetermined  (and  sometimes 
large)  value  of  the  damping  coefficient  6. 


Arguments  must  now  be  made  regarding  contribution  of  the  0(1)  and  (})0 
cross  derivative  terms  in  the  equation  set.  The  goal  is  to  demonstrate 
that  the  eigenvalues  of[L]  ^[W]give  a minimum  At,  so  that  the  0(})  and  (})0 
terms  can  be  dropped  from  consideration.  First,  assume  that  the  contri- 
bution of  these  terms  is  ^ £ -j-  This  effectively  increases  the 
contribution  (results  in  a lower  At)  of  these  terms  since  Ae«A4>.  As  an 
example,  the  matrix[L]*^ [X]is 

[L]"^[X]  = To  0 0 0 oH 
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The  maximum  eigenvalue  of  this  matrix  is  (BX/p)  /(pRe)  and  is  clearly 
less  than  the  eigenvalue  due  to  normal  stress  damping  of  the  [L]'^[W] 
matrix.  Since  the[L]  ^[S]matrix  is  similar,  it  can  be  concluded  that 
sufficient  information  is  known  to  postulate  a complete  stability  criteria. 
By  substituting  the  eigenvalues  into  the  appropriate  stability  criteria 
and  labeling  each  resulting  At  as  to  its  source,  the  following  complete 
stability  criteria  is  postulated: 


At  = m min  (At^^^,  At^^^,  At^.^J 


(25) 


where  m £ 1.0,  = min  of  results  from  Equation  21 


A-t-  - pRe  A9^ 

^^nsd  2 1 (2+^7( 


At.,i,„ pRe  Pr  A0 

vise  = mm  •^—^5 

2 7 


where  min  implies  the  minimum  value  of  the  quantity  found  by  searching 
all  mesh  points  in  the  flow  field.  The  factor  m is  theoretically  limited 
to  a value  of  1.0.  However,  uncertainties  and  approximations  in  the 
stability  analysis  sometimes  result  in  an  actual  maximum  value  of  m 
greater  or  less  than  1.0. 

The  above  procedure  was  tested  at  a = 0°  for  Tracy's  conditions  with 
good  results.  However,  the  real  issue  was  whether  the  criteria  would  be 
adequate  for  high  a with  the  use  of  large  amounts  of  normal  stress 
damping.  Therefore,  the  a = 20"  case  at  Tracy's  conditions  was  repeated 
with  the  above  criteria  for  determining  At  in  continuous  use.  Use  of 
the  criteria  reduced  the  number  of  time  steps  necessary  for  convergence 
to  the  fifth  significant  digit  from  1700  to  960.  The  factor  m for  this 
run  was  1.0  for  the  first  300  time  steps  and  1.2  thereafter.  Instability 
of  the  solution  occurred  for  m = 1.4.  The  values  of  the  normal  stress 
damping  coefficient  (B)  were  comparable  with  those  used  in  the  original 
run.  In  order  to  ensure  that  the  viscous  stability  limit  would  not  be 
violated  when  using  very  high  values  of  normal  stress  damping,  a run  was 
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made  at  Tracy's  conditions  and  a = 0°  with  (5  = 4000.  This  set  of 
conditions  was  unstable  at  a value  of  At  determined  from  Equation  21. 
However  the  use  of  the  above  criteria  reduced  the  value  of  At  sufficiently 
to  allow  the  run  to  be  continued  (although  at  a very  low  At). 

In  summary,  stability  criteria  were  postulated  and  confirmed  which 
reduced  the  number  of  time  steps  required  for  convergence  by  44%  in  an 
example  calculation.  The  procedure  was  also  shown  to  account  for  the 
effects  of  normal  stress  damping  on  stability. 
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